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ABSTRACT 


The bearingless rotorcraft offers reduced weight, less complexity 
and superior flying qualities. fThe bearingless rotors are presently 
being developed and it is most likely that the next generation rotor- 
craft would be equipped with these rotors. Almost all practical 
designs of bearingless rotors include multiple load paths and the 
one that was flight tested by the Boeing Vertol has three load paths. 
The determination of natural vibration characteristics is basic to any 
dynamic design and they form the basis for all practical aeroelastic 
stability analyses of rotor blades. Almost all the current industrial 
structural dynamic programs of conventional rotors which consist of 
single load path rotor blades employ the transfer matrix method 
because this method is ideally suited for one-dimensional chain-like 
structures. In this report, this method is extended to multiple 
load path rotor blades without resorting to an equivalent single 
load path approximation. Unlike the conventional blades, it is 
necessary to introduce the axial-degree-of- freedom into the solution 
process to account for the differential axial displacements in the 
different load paths. With the present extension, the current 
rotor dynamic programs can be modified with relative ease to account 
for the multiple load paths without resorting. to the equivalent 
single load path modeling. The results obtained by the transfer 
matrix raetnod are validated by comparing with the finite-element 
solutions. A differential stiffness matrix due to blade rotation 


is derived to facilitate the finite-element solutions. 
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NOMENCLATURE 


Area of cross-section 
Cross-section integrals, reference 8 
Cross-section integrals, reference 8 

Mass centroid offset from elastic axis, positive when in 
front of elastic axis 

Distance between area centroid and elastic axis, positive 
for centroid forward 

Young's modulus 

Shear modulus 

Cross-section moment of inertia about y-axis 
Reference moment of inertia for non-dimensionalization 
Cross-section moment of inertia about z-axis 
Torsional rigidity 

Polar radius of gyration of cross-sectional area effective 
in carrying tensile stresses about elastic axis 

Mass radius of gyration of blade cross section 



m 




Cross-section integrals, reference 8 


Length of the load paths 
Mass per unit length 

Reference mass for non-dimensionalization 
Twisting moment about x-axis 

Bending moments about y and z axes, respectively 

Number of load paths 

Axial force 

Radius of the rotor 


Tension 



V 


[T] 

[T 1 ] 

U, V, w 

V , V 

y z 

x, y, 2 
{ 2 } 

V 
0 
4 > 

u) 

a 

0 

pc 

Superscripts 


i 


Transfer matrix of the blade 
Transfer matrix of the ith load path 

Elastic displacements in the x, y, 2 directions, respectively 

Shear forces along y and z directions, respectively 

Mutually perpendicular axis system with x along the 
undeformed blade and y towards the leading edge 

State vector 

Slope of deflection curve normal to plane of rotation 
Slope of deflection curve in the plane of rotation 
Pretwist angle 

Elastic twist about the elastic axis 
Frequency of vibration 
Blade rotational speed 
Precone angle 


Differentiation with respect to x 
Differentiation with respect to time 
Non-dimensional quantity 

Quantities corresponding to the ith load path 


Subscripts 

1 Quantities at the root 

2 Quantities at the clevis 

3 Quantities at the blade tip 

i Quantities corresponding to the ith load path 



I. INTRODUCTION 


The bearingless rotorcraft offers simplicity of the design and 
superior flying qualities. The original purpose of introducing hinges 
is to relieve the blades from high stresses and this is an important 
design problem for bearingless rotors and a solution to this problem 
lies in the optimization of blade root stiffness distribution and 
the use of advanced composite materials. The bearingless rotor 
technology was successfully applied by Sikorsky in Blackhawk and 
S-76 helicopters for their tail rotors (Ref. 1). Boeing Vertol 
built the first successful bearingless main rotor (Refs. 2 and 3) 
and this flew first in 1978. During a recent study (Ref. 4) for 
concept definition of the Integrated Technology Rotor/Flight Research 
Rotor (ITR/FRR), thirty- three hub concepts were proposed. Twenty-one 
out of these thirty-three concepts were bearingless designs and hence 
it is very likely that the next generation rotorcraft would be equipped 
with a bearingless rotor. 

Basically, the structural design of bearingless rotors include 
multiple load paths and also some kinematic couplings are introduced 
intentionally through design for specific reasons. It is evident from 
Refs. 4 to 7 that several flexbeam and pitch control configurations 
are possible for bearingless main rotors and the variations in the 
configuration have significant effects on the aeroelastic stability 
of the rotors. Therefore, it is necessary to have accurate analytical 
methods to predict the natural vibration and aeroelastic stability 
characteristics of rotor blades including the multiple load paths. 
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The determination of natural vibration characterisitcs is basic 
to any dynamic design of the rotor system. Also, they form the basis 
for almost all practical aeroelastic stability analyses of rotor blades. 
Almost all the current structural dynamic programs of single load path 
rotor blades employ the transfer matrix method because the method is 
ideally suited for one- dimensional chain-like structures. In this 
report, this method is extended to multiple load path rotor blades with- 
out resorting to the equivalent blade modeling. Any equivalent single 
load path approximation may not simulate completely the dynamic behavior 
of multiple load path systems. With the present extension, the current 
rotor dynamic programs can be modified with relative ease to account 
for the multiple load paths without resorting to the equivalent single 
load path approximation. The results obtained by the transfer matrix 
formulation are validated by comparing with the finite-element solutions. 
The finite-element solutions are generated by adding a differential 
stiffness matrix associated with the rotation to the non-rotating 
stiffness matrix of the blade. To facilitate these calculations, a 
general purpose differential stiffness matrix due to rotation is 
derived analytically. This matrix will be very useful both for the 
generation of finite-element solutions and validation of the 


transfer matrix solutions. 
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II. BASIC EQUATIONS 


The nonlinear equations of motion for the elastic bending and 

torsion of rotor blades are given below from Ref. 8. For algebraic 

* * * 

conciseness, the terms e^, , C^, and are treated 

as zero. However, they do not affect the general nature of the 
formulation presented here. Also, the aerodynamic terms L u> I> v> 
and L w are omitted from the equations of motion. To account for the 
differential axial displacements in the load paths the axial degree- 
of-freedom is included in the equations of motion. 


- {EA(u' + y + y )}’- 2ftmv 

+ mu - ft^u = C^mx 

- (Tv 1 )' + {[El , cos 2 (0+<J>) + EI__, sin 2 (0-Hj>) ] v" 

z y 

+ (El , - El , ) cos(0+<|>) sin(0-H»)w"}" 
z y 

***** * • 

+ 2flmu + mv - me<}> sin0 - 2mefi(v' cos0 + w'sin0) 

2 

- mft [v + e cos(0+<|>)] - 2mQ8p C w 

- {me [ft 2 x cos(8+<j>) + 2ftv cos0]}' = 0 (2) 


- (Tw’) + {(El , - El ,) cos (0+cfr) sin (0+$) v" 

z y 

+ [El , sin 2 (d+<p) + El , cos 2 (0+<*>)] w"}" 
z y 

+ mw + me<b cos0 + 2mfi8 v 

pc 

2 2 

- {me [fi x sin (0+$) + 2S2v sin©]}' = - mfi B x 


( 3 ) 
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— ( G J ' ) * + (El , - El ,)[(w" 2 - v" 2 ) cos0 sin0 
z y 

+ v" w" cos20] + mk 2 '<j> + mfl 2 <() (k 2 - k 2 ) cos20 

m m 2 m^ 

0 •• 9 

+ me [fix (w'cos0 - v'sin0) - (v - fl v) sin0 


•• 2 2 2 2 
+ w cos0] = - 0 m (k -k ) cos0 sin0 - meJl 0 x cos0 

m 2 m^ pc 


(4) 


,2 ,2 

where T =» EA (u* + ^ + j ) 


(5) 


Equations (1) to (4) are coupled nonlinear partial differential equations in 
variables u, v, w and <j>. Substitute Eq. (5) into Eq. (1) and integrate as 
shown below 


T(x) 


T'(x)dx = - f n“mxdx - 


J 

o 


m(2Qv - u + ft u)dx + k^ + k 2 (6) 


where k^ and k 2 are arbitrary constants. 

The boundary condition for tension T(x) is given by 


T(R) = 0 


(7) 


Substitute Eq. (6) into Eq. (7) 


R 


R 


ft mxdx - 


m(2ftv - u + u)dx + k^ + k 2 = 0 


( 8 ) 


To satisfy Eq. (8) , take k^ and k 2 as 
R 


k l = 


2 

(2 mxdx 


(9) 


k 2 


m(2ftv - u + ft u)dx 


( 10 ) 


Substitute Eqs. (9) and (10) into Eq. (6) 
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T(x) 


fl mxdx + 


m(2ftv - u + Si u)dx 


From Eq. (5), u' is given by 


, = T_ v' 2 w' 2 

1 “ EA " 2 "2 


(ID 


( 12 ) 


Substitute Eq. (11) into Eq. (12) 


u'(x) 


EA 


mxdx + — 


1_ 

EA 


R 

2 

m(2ftv - u + SI u)dx 
x 




(13) 


Integrate Eq. (13) with respect to x 

x R x R 


ST f f iff -2 

u(x) * — mxdx + — J J m(2ftv - u + Si u)dx 


o x 


o x 


J. 

2 


(v 


.2 . ,2 

' + w 


) dx 


(14) 


The constant integration is zero because u(0) = 0. Differentiate Eq. (14) 
with respect to time 


x R 

/■ r 

rs 

m(2Siv - u + Si u)dx 

o x 



x 

( v * v * + w'w')dx 

j 

o 


(15) 


Now, by substituting Eqs. (11) and (15) into Eqs. (2) and (3) for the under- 
lined terms, two equations can be obtained in terms of v, w and <f>. These two 
equations plus Eq. (4) are the necessary equations to solve for the unknowns 
v, w and 4>. Once v and w are known, u can be obtained from Eq. (15) and T(x) 
can be obtained from Eq. (11). In essence, the bending and torsional equa- 
tions are decoupled from the axial equation by the above formulation. 

For determination of natural vibration characteristics, one is interested 
in the linear, homogeneous, undamped equations of motion. The process of 
obtaining these equations involves the following steps: 



1 . 


6 . 


Substitution of Eqs. (11) and (15) into Eqs. (2) and (3). 

2. Substitution of the following relations for small values of $. 

cos(9 + $) = cos0 - <J> sin0 

cos 2 (0 + $) = cos 2 0 - $ sin29 

sin(0 + <j>) = sin0 + <f> cos0 

2 

sin (d + <p) = sin20 + ip cos20 

3. Dropping of nonlinear terms, that is, terms containing the products 
of u, v, w, <p and/or their derivatives. 

4. Dropping of damping type terms, that is, terms containing u, v, w, 

$ and their derivatives. 

5. Dropping of nonhomogeneous terms, that is, terms not containing 
u, v, w, <p and/or their derivatives. 

The above steps reduces Eqs. (1) to (4) to the following equations for simple 
harmonic time dependency of u, v, w and <p with frequency uj. 

(EAu')' + (u) 2 + n 2 )mu = 0 (16) 

- (Tv')' + {(El cos 2 0 + El sin 2 0)v" 

z y 

+ (El - El ) cos0 sin0 w"}" - m 2 m v 

z y 

2 2 
+ a) me<J> sin0 - mfl (v-e sin9<{)) 

+ (mefl 2 x sin0$) ' = 0 (17) 

- (Tw 1 )’ + {(El - El ) cos0 sin0 v" 

z y 

+ (El sin 2 0 + El cos 2 e)w"}" - w 2 mw 
z y 

2 2 

- a) me cos0<j> - (meft x cos0cf>) ' = 0 


( 18 ) 



7 . 


- (GJV)' - to 2 mk 2 <j> + 0 2 m (k 2 - k 2 )cos29<{> 


m 


m 2 


m 


1 


2 2 2 
mefl x (cos0w' - sinSv') + me (to + ) sindv 


- to me cosSw = 0 


(19) 


where 


T - fl 


mxdx 


( 20 ) 


The above equations can be reduced to the following twelve non-dimensional 
first-order differential equations: 


du/dx = (El /EAR 2 )N 
y o 


( 21 ) 


dw/ dx = tj) 


dv/ dx » v 


dij;/dx - C._ M + C,, M 
12 z 11 y 


( 22 ) 

(23) 

(24) 


dv/dx 

d^/dx 


dM /dx 
x 


C 00 M + C„, M 
22 z 21 y 


(El /GJ) M 
y o 


_2 _2 -2 

- to me cos0 w + (to + ft ) me sin9 v 


_2 - -2 

+ fl mxe cos9<J; - 0 mxe sin0 v 


(25) 

(26) 


+ (ft 2 m (k 2 - k 2 ) cos20 - to 2 mk 2 }$ 

m 2 m l m 2 


(27) 


dM /dx 
z 


dM /dx 

y 


Tv - ft 2 mex sin0<£ - V 


- Ttj) - ft 2 mex cos0$ + V 


(28) 

(29) 


dV /dx 

y 


dV /dx 
z 


-2 -2 -- -2 -2 

- (to + ) mv + (to + ft ) me sin0$ 

_2 _2 

- to mw - 3" me cos0<J> 


(30) 

(31) 


_2 -2 

= — (to + !] ) mu 


dN/dx 


(32) 
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where 


x = x/R ; u = u/R ; v = v/R ; w = w/R ; ip = ip ; v = v ; 
i - d> ; M = M R/EI ; M R/EI ; M = M R/EI ; 

y v x x y Q y y Q z z y Q 

N = NR 2 /EI ; V = V R 2 /EI ; V = V R/EI ; i = m/m n ; 

y y y y z z y ° 

3 'o 'o 

. r 2 ,2 . 2 r 2 ,2 /t} 2 r 2 , 2 /T1 2 

e = e/R ;k = k lR ; k =k/R;k =k/R; 

m^ m^ mm 

u 2 = u) 2 m R 4 /EI ; Q 2 = fi 2 m R 4 /EI ; C,. - C.. El ; 

o y o y 11 11 y 

3 o 3 o 3 o 

^12 = C 12 EI y 5 ^22 = C 22 EI y 5 E 21 = C 21 EI y 

3 o o o 


El = Reference Stiffness ; m = Reference Mass 

y, ° 


-1 


o 





11 

C 12 


a il 

a 12 

21 

C 22 

S3 

a 21 

a 22 





_ 


a, . = - (El cos 2 0 + El sin 2 8) 
11 y z 


a. = - a 01 = (El - El ) cos0 sin0 
12 21 y z 

2 2 

a = (El sin 0 + El cos 0) 

22 y z 


1 


T = 


ft 2 mxdx 


Equations (21) to (32) can be arranged into a matrix differential equation of 
the following form: 

{z' (x) } = [A(x) ] (z(x) } (33) 

The transfer matrix relating the state vector at any location to the initial 
state vector is defined by 

{ z(x) } = [T(x) ] (5(0)} 


(34) 
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It can be shown that transfer matrix defined in the the above equation is 
governed by the following equation (Ref. 9): 

[T’(x>] = [A(x) ] [T (x) ] (35) 

[T(0)] = identity matrix (36) 


where [A(x) 1 is defined in Eq. (33). 
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III. NATURAL VIBRATION CHARACTERISTICS 


Almost all the practical bearingless rotor designs include multiple 
load paths. For instance, the bearingless main rotor that was flight tested 
by the Boeing Vertol has three load paths, viz., two fiberglass flexbeams 
and a filament-wound graphite torque tube (Ref. 7). The flexbeams and torque 
tube are connected to the blade through a rigid clevis. The idealization of 
a three load path rotor blade is shown in Fig. 1. 


3.1. Equilibrium Across the Clevis : 

Consider the free-body diagram for the clevis as shown in Fig. 2. Let 

(h , h ) be the location of the ith load path with reference to a coordinate 
y i Z i 

system located at the blade (point ‘O'). 

Force equilibrium requires the following relations to be satisfied: 


n 

N = I N 
i=l 1 


(37) 


n 

E V 
i-1 7 i2 


(38) 


n 

V = E V (39) 

Z 2 i=l z i2 


Moment equilibrium requires the following equation to be satisfied. 

n 


(iM + j M + k M ) - 

~ X„ = V- ~ Z„ ' 


E 

i=l 


(iM + M + k M ) 
- x. 0 - y. n - z. 0 

i2 J i2 i2 


- E (jh + k h ) x (i N,, + j V + k V ) 

- y. ~ z. - i2 ~ y. 0 ~ z.' 

1=1 i l ^i2 i2 


The above equation can be written in the component form as shown below: 
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M = E (M -hV +hV) 


"2 i=l "12 z i y i2 J i “±2 




(40) 


M = E M + h N.« 
y 2 1=1 y 12 2 1 12 


□ 

M = E M -h N. 0 
z 2 1-1 2 12 y l 12 


(41) 


(42) 


Equations (37) to (42) can be arranged into a matrix equation as shown below. 

n 


where 


(f } = E [A ] {f } 
1=1 1 


(f 0 ) T - M M M V V N 0 

2 L X 2 Z 2 y 2 y 2 Z 2 U 


{f 


12 


} T - |m M M V N,„ I 

[_ X i2 Z i2 y i2 y i2 2 12 i2 J 


(43) 


and 


[A ± ] = 


1 

0 

0 

0 


0 

1 

0 

0 


0 0 


0 

0 

1 

0 


-h 


0 

0 

1 


y i 

0 

0 

0 


-h 


3.2. Compatibility Across the Clevis 

Consider the plane of the clevis as shown in Fig. 3. Let (x,0,0) t> e 

the location of the blade and (x, h , h ) be the location of the ith load 

y . z . 

ii 

path before deformation. Assume that the plane of the clevis is rigid. 
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that is, rotates about x, y and z axes without any elastic deformation. 
After deformation, the x-coordinates of the blade and the ith load path 
are given by the second line of the following table. 


State 

Axial Coordinate 

Before 

Deformation 

After 

Deformation 

Blade 

ith Load Path 

X 

X + U£ 

X 

x + u 0 - h K - h v„ 

2 z i 2 v 2 


The axial displacement of ith load path is given by subtracting the x-coordinate 
before deformation from the x-coordinate after deformation as shown below. 

u 12 ' <x + “2 - h 2l *2 ' \ V ' x 

Rearrangement of the above equation yields 

"2 “ “12 + h 2 1 *12 + h 7l “ l 2 (44> 

The other compatibility conditions consistent with the rigid clevis are 

“2 ' "12 : v 2 ' v 12 : *2 = *12 : “2 ' v 12 

The above equations together with Eq. (44) can be arranged into a matrix 
equation of the following form 




16 . 



or 

{d 2 } = [B ± ] {d l2 ) 


(45) 


3.3 Frequency Determinant 

Define the following relations between the state vectors from Fig. 1 
and the definition of transfer matrix 


{z„} = [T] {z,} (46) 

Equation (46) relates the state vectors of the blade at Stations 2 (clevis) 
and 3 (tip) . Rewrite this equation into the following partitioned form. 



where d stands for deflections and f stands for forces. 

Extract the following equation for forces from Eq. (47) 


(47) 


{ f 3 } = [T 21 ] (d 2 > + [T 22 ] {f 2 > 


(48) 
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Similarly, the transfer matrix relation for the ith load path can be written 
as (see Fig. 1) 


{z i2 > - [T 1 ] { Zil } 


Rewrite the above equation into partitioned form as 


+ 

i 1 i 

T 1 i t 1 

i 21 l 22 


The displacement and force vectors can be written in terms of {f by virtue 
of boundary condition {d^} = {0} as shown below. 

{f i2 } - [T* 2 ] {f iX } (51) 

{d i2 } = [T 12 ] tf il } (52) 


From Eq. (52) 


{f il } [T 12 ] {d i2 } 


Substitute Eq. (45) into Eq. (53) 

i - 1 

{f il } = CT 12 ] [B i ] {d 2 } 
Substitute Eq. (54) into Eq. (51) 


. . -1 -1 
{f i2 } = [T 22 3 tT 12 ] [B i ] {d 2 } 


Substitute Eq. (55) into Eq. (43) 


n . -1 -1 

{f 2 > = l [A.] [T 22 ] [Tj 2 ] [B.] (d 2 ) 


Substitute Eq. (56) into Eq. (48) 
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{f 3 > = ( [t 21 ] + 


[ t 22 ] 


n 

E 

1=1 


[A ± ] [Tj 2 ] 




-1 


-1 


[B ± 3 


) {d 2 } 


(57) 


The vector {f^} = {0} by virtue of the boundary conditions and for nontrivial 
solutions of (d 2 ), the determinant of the coefficient matrix should be zero 
and this condition yields the following frequency equation to determine the 
natural frequencies. 


det ([T 21 ] 


[t 22 ] 


n 

E 

i=l 


-1 


[A ± ] [T 22 ] [tJ 2 ] [B i ] 


-1 


(58) 


3.4 Mode Shapes 

The state vector at the clevis from Eq. (47) can be written as 

T 
T 


where 


'll 

1 f 
1 l 12 
i 

21 

“ “T 

I m 

1 T 22 


(59) 


l 

T 1 

1 12 ! 
1 

T 12 


” T 11 

1 

1 

i 

— 

T 

1 12 

i— ^ 

CNI 

IH 

1 

t 22 


T 21 

'7 ” 
1 

1 

1 

H i 
N> 1 


-1 


From the above equation, the displacement vectors at the clevis and the tip 
of the blade can be related as shown below by virtue of the boundary condition 


(f 3 ) « (0}. 


{d 2 } = {d 3 } 


(60) 


Substitute Eq. (60) into Eq. (57) 


where 


{f 3 > = [C] (d 3 ) 


n -1 

[c] = ([t 21 ] + [t 22 ] z [a.] [t 22 ] [ t \ 2 ] 


( 61 ) 


[B.]) [T u ] 
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Assume w^ = 1 arbitrarily and rewrite rows 


2 to 6 of Eq. (61) as 



(62) 


By solving the above equation, u^» v^, an( * ^3 are known and together 

with w^ = 1 the entire (d^) is known. Once, {d^} is known, the state vector 
at the clevis can be determined from Eq. (59) as shown below 



(63) 


Once the state vector at the clevis is known, the deflection vectors in the 
blade and the load paths can be obtained as follows: 

Blade: By definition of the transfer matrix the state vector at any 

location x, is given by 




From the above equation 


{d x } = [T 11 ] {d 2 } + [T 12 ] {f 2 } 

x x " 


(64) 


(65) 


Load path: By definition of the transfer matrix the state vector at 


any location, x, in the load path is given by 



20 . 




From the above equation 


<d > - Ct 12 ] tf } 

X * 


by virtue of the boundary condition ({d^} ° {0})* 


Substitute Eq. (54) into Eq. (67) 

{d lx> - ^ ‘V’ 1 < d 2> 

x 


The mode shapes can be computed from Eqs. (65) and (68). 


(65) 


(67) 


( 68 ) 


3.5 Tension Coefficients 

The preceding formulation assumes that the transfer matrices for the 
blade and the load paths are either known or can be calculated. The trans- 
fer matrices can be calculated if all the coefficients of the differential 
equations of motion are known and tension appears as coefficient in these 
equations. The calculation of tensions in the blade is straight forward 
and it is calculated from Eq. (20). 

The following procedure is employed for calculation of the tension 
coefficients in the load paths. Draw a free-body diagram for the clevis 
as shown in Fig. 4. The tension 'T' applied by the blade to the clevis 
is known from Eq. (20). From the free-body diagram the following equilib- 
rium equation can be obtained 

n 

{b} = Z [A.] (f } 
i=l 1 


( 69 ) 





22 . 


where 


{ b} T = |_0 0 0 0 0 T_J 

t ] is defined in Eq. (43) 


Let [T ] be the transfer matrix of the ith load path corresponding to the 
static case (oi = 0) . By definition, one can write that 


mi 

1 

ml " 

T 

L n 

1 

1 

T 

12 

mi 

T 

1 

mi 

T 

21 

I 

i 

T 

22 


From the above equation 


{f i2 } = [T 22 ] {f il } 


Substitute Eq. (71) into Eq. (69) 


{b> = E [Aj [T^] {f tl } 


From Eq. (70) 


{d i2 } = [Tj 2 ] {f iL } 


by virtue of the boundary condition ({d^} = ^}) 

Substitute Eq. (73) into Eq. (72) 

n . -1 

(b) = Z [A.] tT 22 ] [T n ] {d. 2 } 


Substitute Eq. (45) into Eq. (74) 

n . -1 -1 

{b> = E [A ± ] [T 22 ] [T^ 2 ] [B ± ] {d 2 > 

i=l 


( 75 ) 
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From the above equation 


-1 

{d 2 > = [D] {b} 


(76) 


where 


n i ' i - 1 - 1 

[D] = ^ [A 1 ] [T^] [Tj 2 ] [B i ] 


From Eq. (45) 


-1 

(d i2 } = [B i ] {d 2 > 


(77) 


From Eq. (73) 


( f il > = ^ T 12^ {d i2 } 


(78) 


Substituting Eqs. (76), (77) and (78) into Eq. (71) 

i i - 1 - 1 - 1 

{f i2 > - [T 22 ] [T\ 2 ] [B i ] [D] (b) 


(79) 


The last or sixth row in the above equation is N i2> Once N i2 is known, the 
tension in the ith load path can be obtained from the following equation 


T.(x) 




l 

mx dx + N^ 2 
x 


(80) 


where 


i = span of the load path. 


The above formulation assumes that the transfer matrices (static) of 
the load paths are known which depend on the tensions to be calculated. 
The following iterative scheme is used to solve the problem. The initial 
tensions in the load paths are calculated by assuming that the load paths 
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are coincident with the blade, i.e., h = h =0. The tensions corres- 

y i Z i 

ponding to this case are calculated as follows. The transfer matrix for 
axial motion of a beam for the static case is given by 

1 a 

0 1 


[T(x) ] 


(81) 


where 


dx 

EA 


By definition of the transfer matrix 


u 


i2 


N 


i2 


“ii - 0 


N 


il 


(82) 


From the above equation 


"12 ' a i N il 


N i2 ■ N U 


Consider the following two cases! 


Case 1 : 2 Load Paths 


(83) 


Expanding Eq. (83) 


U 12 = a l N 12 


u 22 a 2 N 22 


( 84 ) 


For coincident nodes (u^ = u 22 ) 


Eq. (84) becomes 
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a l N 12 " a 2 N 22 " 0 


For equilibrium 


N 12 + N 22 = T 


By solving Eqs. (85) and (86) 


N 


12 a 1 + a 2 


N 


22 a^ + a 2 


Case 2: 3 Load Paths 


Expanding Eq. (83) for coincident nodes 


a l N 12 * a 2 N 22 a 3 N 32 


For equilibrium 


N 12 + N 22 + N 32 = T 


By solving Eqs. (88) and (89) 


N 


a 2 a 3 T 


12 


N 


a 3 a l T 


22 


N 


a l a 2 T 


32 


where 


D = (a^ a 2 + a 2 a^ + a^ a^) 


(85) 


( 86 ) 


(87) 


( 88 ) 


(89) 


( 90 ) 
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Now Che above result can be generalized to the n load path case as shown 
below. 


where 



n,i 



n n,k 


£ ( it a.) 

k-1 j=l J 


n,l 

it a, = a. a 9 
j = l J 



(91) 


Now, the tension in the ith load path corresponding to the coincident load 
path case is given by 

a 

c 2 

T^(x) = Q mx dx + N 12 (92) 

o 

The tensions obtained from coincident nodes are used to obtain the trans- 
fer matrices of the noncoincident problem. At the most, one more itera- 
tion may be required to obtain the convergence for practical bearingless 
rotor designs. 
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IV. RESULTS 

4.1 Computer Program 

A general purpose computer program has been developed to determine 
the natural vibration characteristics of rotating multiple load path 
rotor blades based on the transfer matrix formulation. The listing 
of the computer program and the user's instructions are given in 
appendices A and B respectively. The significant features of the 
computer program are as follows. 

1. The non-uniform properties including the pretwist for the 
blades as well as for the load paths can be handled by 
the program. 

2. The coupled flapwise bending, chordwise bending, torsion 
and axial stretching degrees-of-freedom are included in 
the program. In the case of single load path blade the 
axial stretching is decoupled similar to the conventional 
analysis. 

3. No equivalent single load path approximation is made in 
the program. 

4. Maximum number of load paths allowed in the program is three 
but can be extended very easily to "n" number of load paths. 

5. Continuous system model is used in the program and a fourth- 
order Runge-Kutta integration scheme is used to determine 
the transfer matrices. If a discrete model is used just 
like in the case of Bell Helicopter's C-81 program, the cor- 
responding transfer matrices can be used in place of the 
continuous system transfer matrices and the formulation 

is independent of the model. 
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6. The natural frequencies are computed by a frequency scanning 

technique. Sign changes in the values of frequency determinant 
are detected by starting from an initial value and incrementing 
at steps of specified "h" until the required sign changes 
are detected or a final prescribed value is reached. If 
two frequencies are closer than the increment "h", then 
there is a chance of missing those frequencies. In such 
cases, the frequencies are detected by the fact that if 
any three consecutive frequency determinants have the 
same sign and the absolute value of the middle determinant 
is the smallest of the three, then there are two frequencies 
in that range. In this case smaller increments are taken 
to bracket the roots. 

4.2 Numerical Results 


The following data is used for numerical calculations to validate 
the formulation and the computer program. 

1. Radius of the rotor = 260 in. 


2. Distance of the clevis from the root 


= 52 in. 


3. Length of the blade 

4. Rotational speed (ft) 

5. Flapwise bending stiffness (EI^) 

6. Chordwise bending stiffness (El ) 

z 

7. Torsional stiffness (GJ) 

8. Axial stiffness (EA) 

9. Mass per unit length 

10. Built-in- twist 

11. Collective pitch 


= 208 in. 

= 360 RPM 

= 0.2977 x 10 8 lb/in 2 
= 10 x 10 8 lb/in 2 

= 0.2 x 10 8 lb/in 2 

11 2 
- 10 lb/ in 

= 0.0015 lb-sec 2 /in 2 

= 0 deg. 

= 15.026 deg. 
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12. e 

13. mk 2 

m. 


14. mk 2 
m 2 

For two load path blade 



0.0; h = -1.0 in; h 

Z 1 y 2 



3.0 in. 


= -0.6 in. 

= 0.89545 x 10~ 3 lb-sec 2 
=0.04 lb-sec 2 


For three load path blade 

h = 1.0 in; h = 3.0 in; h = -1.0 in; h = -1.0 in; 
y l Z 1 y 2 Z 2 

h = - 2.0 in. 



2.0 in; 


The natural frequencies obtained from the computer program are presented 
in Table 1 for load paths 1 and 2 respectively corresponding to the 
rotational speed ft = 360 RPM. The non-rotating and rotating (ft = 360 RPM) 
natural frequencies for a three load path case are presented in Table 2 
and in this case out-of-plane locations are assumed for the load paths 
for validation of a general case. A mode shape corresponding to the 
two load path case is presented in Table 3. All the results obtained 
from the transfer matrix formulation are compared with the finite- 
element solutions for validation of the present formulation. The 
finite-element solutions are obtained from the NASTRAN program using 
the 80 beam elements. An external stiffness matrix associated with 
the centrifugal forces was computed separately and added to the non- 
rotating stiffness matrix calculated by the NASTRAN program. The 
derivation of the differential stiffness matrix associated with the 
blade rotation is given in the next section. From the results pre- 
sented in Tables 1 to 3, it is clear- that the trnasfer matrix formulation 
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Table 1. Comparison of natural frequencies, single and two load 
path cases, rad/sec 



Single Load Path Case 

Two Load Path Case 

Mode 

No. 

Transfer 
matrix method 

Finite-element 

method 

Transfer 
matrix method 

Finite-element 

method 

1 

36.7738 F 

36.8062 F 

43.1129 F 

43.1412 F 

2 

48.1092 C 

48.1511 C 

56.5803 C 

56.6239 C 

3 

104 .9309 F 

104.9170 F 

122.1807 F 

122.1775 F 

4 

138.2931 T 

138.2254 T 

152.4027 T 

152.3739 F,T 

5 

202.4001 F 

202.3235 F 

226.2707 F 

226.2831 F 

6 

280.5927 C 

280.4453 C 

307.8609 C 

307.7620 C 

7 

336.3352 F 

336.1560 F 

330.6756 F 

330.7502 F 

8 

402.5505 T 

402.5403 T 

431.3825 T 

431.3955 T 

9 

507.5868 F 

507.2334 F 

506.2530 F 

506.3303 F 

10 

669.1642 T 

667.9858 F,T 

669.4276 T 

668.4175 T 


F = Predominantly flapwise bending 
C = Predominantly chordwise bending 
T = Predominantly torsion 
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Table 2. Comparison of natural frequencies, non-coplanar three load 
paths case 


fl = 0 0 = 360 RPM 

Mode Transfer Finite-element Transfer Finite-element 

No. matrix method method matrix method method 


11.4292 F 


11. 4292 F 


44. 8430 F 


44.8613 F 


66.1632 C 


66.2114 C 


72. 8896 C 


73.0002 C 


70.1983 F 


70.2009 F 


124. 3887 F 


124.4651 F 


153. 6332 T 


153.6552 T 


158. 0232 T 


157.9025 T 


18.1. 4035 F 


181.3978 F 


23 2. 6334 F 


232.6989 F 


264. 3820 F 


264.4106 F 


326.3415 F 


326.3315 F 


405.9868 C 


406.1411 C 


418.1878 C 


418.1265 C 


418.2305 F 


418.1628 F 


447.2150 T 


447.6895 T 


445.6527 T 


445.9648 T 


498.3820 C 


498.3477 C 


665.0317 F 


669.7281 T 


668.9731 T 


F = Predominantly flapwise bending 
C = Predominantly chordwise bending 


T = Predominantly torsion 












Table 3. a. Comparison of mode shapes, two load path 
case, third mode. Omega ~ 360 rpm. 

(Flapwise Deflection) 


1 

1 

1 

x/R 

1 

1 

Transfer Matrix 

Method 

[ Finite Element 

j Method 

1 

1 

i 

1 

1 


1 

1 

Load Path 
I 

| Load Path 

1 II 

1 

| Load Path 
1 I 

| Load Path 
1 II 

1 

1 

1 

1 

1 

0.00 

1 

1 

0.0 

| 0.0 

| 0.0 

| 0.0 

1 

1 

l 

1 

0.10 

1 

-0.0743 

1 

| -0,0744 

| -0.0742 

1 

| -0.0743 

1 

1 

1 

1 

1 

1 

0.160 

1 

1 

1 

-0.1315 

i 

| -0.1318 

1 

| -0.1311 

i 

| -0.1315 

1 

1 

1 

1 

1 

1 

1 

0.200 

1 

1 

-0 . 1468 

1 -0.1466 

1 

1 

I 

1 

1 

1 

0.312 

1 

1 

1 

-0.2839 

j -0.2835 

1 

1 

1 

1 

1 

1 

0.408 

1 

1 

1 

-0.4600 

| -0.4591 

1 

1 

1 

1 

1 

1 

0.504 

1 

1 

1 

-0.5732 

| -0.5720 

1 

1 

1 

1 

1 

1 

0.600 

1 

1 

1 

-0.5639 

| -0.5628 

1 

1 

1 

1 

1 

0.712 

1 

1 

1 

-0.3533 

S -0.3521 

1 

1 

1 

1 

1 

1 

0.808 

1 

1 

1 

+0.0008 

| +0.0019 

1 

1 

1 

1 

1 

1 

0.904 

1 

1 

1 

+0.4747 

| +0.4754 

1 

1 

1 

1 

1 

1.00 

1 

1 

1 

+1.0000 

| +1.0000 

1 

1 

1 



x/R 


0.00 

0.10 

0.160 

0.200 

0.312 

0.408 

0.504 

0.600 

0.712 

0.808 

0.904 

1.00 


Table 3. b. Comparison of mode shapes, two load path 
case, third mode, Omega = 360 rpm. 

(Chordwise Deflection) 


Transfer Matrix 

1 

Finite 

Element 


Method 


1 

Method 


Load Path | Load Path 

1 

Load Path 

| Load Path 


I | 

II 

1 

I 

1 II 

1 


0.0 | 
1 

0.0 

1 

0.0 

| 0.0 
1 


1 

0.0210 | 
i 

0.0209 

1 

1 

0.0210 

| 0.0210 
1 


1 

0.0380 | 

1 

0.0379 

1 

1 

1 

0.0380 

i 

| 0.0381 

1 


0.0436 


1 

1 

0.0437 


0.0835 


1 

1 

1 

0.0838 


0.1308 


1 

1 

1 

0.1313 


0.1590 


1 

1 

0.1598 


0.1524 


1 

1 

I 

0.1538 


0.0890 


1 

1 

i 

0.0910 


0.0132 


1 

1 

1 

-0.0019 


-0.1482 


1 

1 

I 

-0.1446 


-0.2972 


1 

1 

-0.2927 




Table 3. c. Comparison of mode shapes, two load path 
case, third mode. Omega = 360 rpm. 

(Torsional Deflection) 


Transfer Matrix 
Method 


x/R 


Load Path 
I 


Load Path 
II 


Finite Element 

Method 


Load Path 
I 


Load Path 
II 


0.00 

i o.o i 

S 1 

0.0 

1 

1 

0.0 | 

0.0 

0.10 

i i 

| 0.0015 | 

0.0015 

1 

1 

I 

1 

0.0016 | 
1 

0.0016 

0.160 

1 1 

| 0.0023 | 

1 1 

0.0023 

1 

1 

1 

0.0026 | 
1 

0.0026 

0.200 

1 

J 0.0029 

1 


1 

1 

0.0032 


0.312 

1 

| 0.0061 
I 


1 

1 

i 

0.0064 


0.408 

1 

| 0.0086 
I 


1 

1 

I 

0.0089 


0.504 

1 

| 0.0109 


1 

1 

I 

0.0122 


0.600 

1 

| 0.0126 
1 


1 

I 

i 

0.0130 


0.712 

1 

| 0.0140 

I 


1 

! 

i 

0.0145 


0.808 

1 

| 0.0147 

i 


1 

0.0153 


0.904 

[ 0.0150 

1 


1 

f 

i 

0.0156 


1.00 

i 

| 0.0151 


l 

I 

0.0157 
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yields very accurate results for multiple load path rotor blades. In fact, 
the slight discrepancies noticed in these results is due to the approxima- 
tions made in the finite-element solutions, for instance, a constant tension 
is assumed within each finite-element and actual variation is used in the 
transfer matrix solution. 


4.3 Differential Stiffness Matrix Due to Rotation 

The work done by centrifugal forces due to combined flapwise bending, 
chordwise bending, torsion and axial stretching is given by (Ref. 10) 

i l 

P-i J T(v ' 2 + w’ 2 )dx - j ft 2 (u 2 + v 2 )dx 


1 

2 


ft me (cosg - $sinB)v dx 


ft m {xe (v' cosg + w* sing) 


+ <J>xe (— v 1 sing + w f cosg)}dx 


4 


ft 2 m (k 2 - k 2 ) (sin2g + <pcos2B) 4> dx 

m2 m^ 


(93) 


Consider the following displacement functions for flapwise bending, chord- 
wise bending, torsion and axial displacement in terms of the shape functions 


w(x) = w. (1 - 3x 2 /£ 2 + 2 y? / 9 ?) + 0 (-x + 2x 2 / 1 - v? / £ 2 ) 

1 y l 

2 9 2 12 

+ w_ (3x / 2, -2 x /£ ) + 0 (x/Jl-x/i) 

2 y 2 


( 94 ) 
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v(x) = v. (1 - 3x 2 /£ 2 + 2x 3 /2) + 0 
X z 

+ V- (3x 2 /2 2 - 2x 3 /2 3 ) + 9, 

2 Z 2 

<j»(x) = 4> 1 (1 - x/2) + <)> 2 (x/ 1 ) 

u(x) = u^ (1 - x/£) + u 2 (x/2) 

Substitute Eqs. (94) to (97) into Eq. (93) and 
tion into the following form 


1 


(x - 2x 2 /£ + x 3 / £ 2 ) 

(-x 2 /£ + x 3 /£ 2 ) (95) 

(96) 

(97) 

arrange the resulting equa- 


P = \ (d} T [k] {d} 


(98) 


where 



u i v i ”i +i 




U 2 V 2 W 2 


0 0 



The matrix [k] in Eq. (98) is the differential stiffness matrix due to rota- 
tion. The final form of this matrix is given in Appendix C. 
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APPENDIX A 



C irirJ e k&irfrii Iviririri. it hit ic!ci r !rkicirlcirlr!r!ririrlrlrk - !r&cirir!ri c !rk^r!r!r!rfri r !riririciririr ir ! r ! r k 4r ! rir!ri; 

C THIS PROGRAM COMPUTES THE NATURAL VIBRATION CHARACTERISTICS OF 
C MULTIPLE LOAD PATH ROTOR BLADES (UPTO THREE LOAD PATHS) 

C UNDERGOING BENDING - BENDING - TORSION - AXIAL VIBRATIONS. 

C 

C 

C DECLARATION STATEMENTS 

C - 

C 

IMPLICIT REAL*8 (A - H, 0 - Z) 

REAL*8 MASS,MASS1(3,21) ,MASS2(101) , KM1S1(3,21) 

REAL*8 KM2S1(3 ,21) , KM1S2(101), KM2S2(101) 

DIMENSION STA1(21), STA2(101), EIY1(3,21), 

1 EIY2(101) , EA1(3 ,21) , EA2(101), EIZ1(3,21), 

1 EIZ2(101), GJ1(3 ,21) , GJ2 (101) , El(3,21), 

1 E2(101) , BETA1(3 , 21) , BETA2(101), tfl(3,ll), 

1 W2(51), Vl(3,ll) , V2(S1), PHI1(3,11) , PHI2(51), 

1 SLl(ll) , SL2(51) , STA(lOl) , FD(3,2), 

1 U(3) , BB(12) , AT ( 3 ) , TEN1(3,21), 

1 TEN2(101) , Dll(3,21) , D12(3,21), 

1 D13(3,21) , D14(3 ,21) , D1S(3,21), 

1 D16 (3 ,21) , D17 (3,21) , D18(3,21), 

1 D19(3,21) , D 191(3 ,21) , D192(3,21), 

1 D21(101), D22(101) , 023(101), 

1 D24(101), D25 (101) , 026(101), 

1 D27 (101) , D28 (101) , D29(101), 

1 D291(101), D292(101) , FREQEN(IO) , 

1 TF1(3,12,12) , TF2(12 , 12) , XA(6,6), XB(6,6), XC(6,6), 

1 XD(6 ,6) , XE(6,6) , XR(6,6) 

COMMON/X1/FREQEN , H, HI, H2, IJK 
C0MM0N/X2 /FACT 
C0MM0N/X3/STA, NS 
C0MM0N/X4/NPATH , ICPL 

COMMON/X5/EIY1 ,TEN2 ,EA1 ,EA2 ,D11 ,D21 ,D12 ,D22 ,D13 , 

1 D23,D14,D24,D15 ,D25 ,D16 ,D26 ,D17 ,D27 ,D18,D28,D19 , 

1 D29,D191,D291,D19Z,D292,0MEGAN 

C0MM0N/X6/SL1 , SL2 , CPM , FRE .HERTZ .BLANK , DOT , STAR , J1 , IPLOT 
COMMON/X7/ HH1, HH2 
C0MM0N/X8/ BB, IND 
COMMON/X9/ FD 
COMMON/XIO/ TF1 , TF2 
C0MM0N/X11/ XA, XB, XC, XD, XE , XR 
COMMON/X12/CON 
C 


A- 1 


noon 


THIS SECTION READS THE DATA OF THE SYSTEM 


READ(20, 100) NPATH, ISTAGE, IPLOT, NS1, NS2 
READ (20, 105) SPAN1, SPAN2, SCH, OMEGA 
READ(20,105) (STAl(I), 1=1, NS1) 

READ (20, 105) (STA2(I), 1=1, NS2) 

DO 400 1=1, NPATH 

READ(20 , 105) (MASSl(I.J), J=1,NS1) 

READ (20, 105) (EIY1(I,J), J=1,NS1) 

READ(20,105) (EIZ1(I,J), J=1,NS1) 

READ(20, 105) (GJ1(I,J), J=1,NS1) 

READ(20 , 105) (E1(I,J), J=1,NS1) 

READ (20, 105) (BETAl(I.J), J=1,NS1) 

READ(20,105) (KM1S1(I,J), J=1,NS1) 

READ (20, 105) (KM2S1(I,J), J=1,NS1) 

READ(20,105) (EA1(I,J), J=1,NS1) 

400 CONTINUE 

READ (20, 105) (MASS2(I), 1=1, NS2) 

READ (20 , 105) (EIY2(I), 1=1, NS2) 

READ (20 , 105) (EIZ2(I), 1=1, NS2) 

READ (20 , 105) (GJ2(I), 1=1, NS2) 

READ (20, 105) (E2(I), 1=1, NS2) 

READ(20 , 105) (BETA2(I) , 1=1, NS2) 

READ (20 , 105) (KM1S2(I), 1=1, NS2) 

READ(20,105) (KM2S2(I), I=1,NS2) 

READ (20 , 105) (EA2(I), 1=1, NS2) 

READ(20 , 105) HI, H, H2 

IF (NPATH. GT.l) READ(20,105) ((FD(I,J),J = 1,2), I = 1, NPATH) 
IF (NPATH. GT.l) READ (20, 100) ITR 
IF ( ISTAGE. NE.l) READ(20,100) NF 
IFUSTAGE.EQ.4) READ(20,105) (FREQEN(J) , J = 1,NF) 
IF(IPLOT.EQ.l) READ (20 , 110) BLANK, DOT, STAR 


A-2 



o o o o 


THIS SECTION PRINTS THE DATA OF THE SYSTEM 


SPAN = SPAN1 + SPAN 2 
WRITE (22 ,200) 

WRITE (22, 205) 

WRITE (22, 347) ISTAGE 

WRITE (22, 348) IPLOT 

IF (ISTAGE .NE. 1) WRITE (22, 210) NF 

WRITE(22,310) HI 

WRITE(22,215) H 

WRITE(22,220) H2 

WRITE (22, 225) SPAN 

WRITE(22,230) SPAN1 

WRITE (22, 235) SPAN2 

WRITE (22, 236) SCH 

WRITE (22, 240) OMEGA 

WRITE (22, 245) NS1 

WRITE(22 ,250) NS2 

WRITE (22, 3 15) NPATH 

DO 405 1=1, NPATH 

IF (I .EQ. 1) WRITE (22, 255) 

IF (I .EQ. 2) WRITE(22,260) 

IF (I .EQ. 3) WRITE (22, 265) 

WRITE (22, 270) 

WRITE(22 ,275) (STA1(J), J=1,NS1) 

WRITE (22, 280) 

WRITE (22 ,275) (MASSl(I.J), J=1,NS1) 

WRITE (22, 285) 

WRITE(22 ,275) (EIY1(I,J), J=1,NS1) 

WRITE (22, 350) 

WRITE (22, 275) (EIZ1(I,J), J=1,NS1) 

WRITE (22, 35 2) 

WRITE(22,275) (GJ1(I,J), J=1,NS1) 

WRITE (22, 354) 

WRITE (22 ,275) (E1(I,J), J=1,NS1) 

WRITE (22, 35 6) 

WRITE(22 ,275) (BETA1(I,J), J=1,NS1) 

WRITE (22, 358) 

WRITE(22 ,275) (KM1S1(I,J), J=1,NS1) 

WRITE (22, 360) 

WRITE(22,275) (KM2S1(I,J), J=1,NS1) 

WRITE(22,286) 

WRITE(22,275) (EA1(I,J), J=1,NS1) 

IF (NPATH .GT. 1) WRITE(22 ,321) I, FD(I,1), FD(I,2) 
IF (NPATH .GT. 1) WRITE(22,375) ITR 
405 CONTINUE 

WRITE (22, 290) 

WRITE (22 ,270) 

WRITE (22 ,275) (STA2(I), 1=1, NS2) 

WRITE (22 ,280) 

WRITE (22 ,275 ) (MASS2(I), 1=1, NS2) 

WRITE (22, 285) 

WRITE(22,275) (EIY2(I), 1=1, NS2) 

WRITE (22, 350) 
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WRITE (22, 275) 
WRITE (22, 352) 
WRITE (22, 27 5) 
WRITE (22, 354) 
WRITE (22, 275) 
WRITE (22, 35 6) 
WRITE (22, 275) 
WRITE (22, 35 8) 
WRITE (22, 275) 
WRITE (22, 360) 
WRITE (22, 275) 
WRITE (22, 343) 
WRITE(22,275) 
WRITE (22, 205) 


(EIZ2(I), I 
(GJ2(I) , I 
(E2(I) , I 
(BETA2(I) , I 
(KH1S2(I) , I 
(KM2S2(I), I 
(EA2(I) , I 


=1,NS2) 
=1,NS2) 
=1,NS2) 
=1 ,NS2) 
=1,NS2) 
=1,NS2) 
=1 ,NS2) 



noon 


THIS SECTION INTERPOLATES THE DATA 


WRITE (22,200) 

WRITE (22, 205) 

WRITE (22, 295) 

NS = NS1 

DO 411 J = 1, NS1 
STA(J) = STA1(J) 
411 CONTINUE 


HH = SPAN1/20 . 0 


420 I = 1 

, NPATH 

DO 410 J 

= 1, NS1 

D21(J) 

= MASSl(I.J) 

D22( J) 

= EIY1(I , J) 

D23(J) 

* EIZ1(I , J) 

D24(J) 

= GJ1(I,J) 

D25(J) 

= El(I.J) 

D26(J) 

a BETAl(I.J) 

D27 ( J) 

a KM1S1(I , J) 

D28(J) 

a KM2S1(I, J) 

TEN2(J) 

= EA1(I,J) 


410 CONTINUE 

CALL INTPOL(21 , D21, HH) 
CALL INTP0L(21 , D22, HH) 
CALL INTPOL(21 , D23 , HH) 
CALL INTP0L(21 , D24, HH) 
CALL INTP0L(21 , D25 , HH) 
CALL INTPOL(21, D26 , HH) 
CALL INTP0L(21, D27 , HH) 
CALL INTPOL(21 , D28, HH) 
CALL INTPOL ( 2 1 , TEN2 , HH ) 

DO 415 J = 1, 21 
HASSl(I.J) = D21(J) 
EIYl(I.J) = D22 ( J) 
EIZl(I.J) = D23 ( J) 
GJl(I.J) = D24( J) 

E1(I, J) = D25 (J) 
BETAl(I.J) = D26(J) 

KM1S1 (I , J) = D27 (J) 
KM2S1(I , J) = D28 ( J) 
EAl(I.J) = TEN2(J) 

415 CONTINUE 

IF (I .EQ. 1) WRITE (22, 255) 
IF (I .EQ. 2) WRITE(22,260) 
IF (I .EQ. 3) WRITE (22 ,265 ) 
WRITE (22, 280) 


WRITE (22, 275) 
WRITE (22, 285) 

(D21( J) , 

J=1 ,21) 

WRITE(22,275) 
WRITE (22 , 350) 

(D22( J) , 

J= 1 , 2 1 ) 

WRITE (22, 275) 
WRITE(22,352) 

(D23 ( J) , 

J=l,21) 

WRITE(22,275) 
WRITE (22, 354) 

(D24( J) , 

J=l,21) 
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WRITE(22,275) (D25(J), J=l,21) 

WRITE(22,356) 

WRITE(22,275) (D26(J), J=l,21) 

WRITE (22, 358) 

WRITE(22,275) (D27(J), J=l,21) 

WRITE (22, 360) 

WRITE(22,275) (D28(J), J=l,21) 

WRITE (22, 286) 

WRITE(22,275) (TEN2(J), J=l,21) 

420 CONTINUE 
NS = NS2 

DO 421 J = 1, NS2 
STA(J) = STA2(J) 

421 CONTINUE 

HH = SPAN2/100 .0 

CALL INTPOL(101, MASS2, HH) 

CALL INTP0L(101, EIY2, HH) 

CALL INTPOL(101, EIZ2, HH) 

CALL INTP0L(101, GJ2, HH) 

CALL INTPOL(101, E2, HH) 

CALL INTP0L(101, BETA2 , HH) 

CALL INTP0L(101, KH1S2, HH) 

CALL INTP0L(101, KM2S2, HH) 

CALL INTPOL(101, EA2, HH) 
WRITE(22 ,300) 

WRITE (22, 280) 

WRITE(22 ,275) (MASS2(J), J=l,101) 
WRITE (22, 285) 

WRITE(22,275) (EIY2(J), J=l,101) 

WRITE (22, 350) 

WRITE (22,275) (EIZ2(J), J=l,101) 

WRITE (22, 352) 

WRITE(22,275) (GJ2(J), J=l,101) 

WRITE (22, 354) 

WRITE(22,275) (E2(J), J=l,101) 

WRITE (22, 35 6) 

WRITE(22 ,275) (BETA2(J) , J=l,101) 
WRITE (22, 358) 

WRITE (22,275) (KM1S2(J), J=l,101) 
WRITE (22, 360) 

WRITE (22 ,275) (KM2S2(J), J=l,101) 
WRITE (22, 343) 

WRITE (22,275) (EA2(J), J=l,101) 


A-6 



n n n n o 


/ 


THIS SECTION NONDIMENS IONALIZES THE DATA AND DETERMINES THE 
COEFFICIENTS OF THE FIRST ORDER DIFFERENTIAL EQUATIONS 


422 


425 


426 


PI = 4.0 * DATAN ( 1 . OD+OO ) 

OMEGA = OMEGA * PI/30.0 

OMEGAS = OMEGA * OMEGA 

IF (NPATH .EQ. 1) FD(1,1) = 0.0 
IF (NPATH .EQ. 1) FD(1,2) = 0.0 
ICPL = 0.0 
DO 422 I = 1, NPATH 

DO 422 J « 1,2 

IF(FD(I,J) .GE. 1.0E-10) ICPL = 1 
CONTINUE 

STEP = SPAN2/ 100.0 
X = SPAN1 
H5 = STEP/24.0 
DO 425 1=1, 101 

D21(I) = MASS2 (I) * H5 * X 
X = X + STEP 

CONTINUE 
TEN2C101) = 0.0 

TEN2(100) = ( 9.0 * D21(101) + 19.0 * 

1 D21(100) - 5.0 * D21(99) + D21(98) ) * OMEGAS 
DO 426 I =2,99 

J =101-1 

TEN2(J) = TEN2(J+1) + ( -D21(J+2) + 13.0 * 

1 ( D2KJ+1) + D21(J) ) - D21(J-1) ) * OMEGAS 
CONTINUE 

TEN2(1) = TEN2 (2) + ( D21(4) - 5.0 * 

1 D21 (3) + 19.0 * D21(2) + 9.0 * D21(l) ) * OMEGAS 
EIY = EIY2(1) 

MASS = MASS2(1) 

SPANS = SPAN * SPAN 

FA = SPANS /EIY 

FACT = DSQRT (FA * SPANS * MASS) 

CON = SPAN / SCH 

OMEGAN = OMEGAS * FACT * FACT 

HH = SPAN1 / (20.0 * SPAN) 

DO 4251 I = 1, NPATH 


X 

DO 4251 J 
KM1S1(I , J) 
KM2S1(I , J) 
EA1(I, J) 
PITCH 
CO 
SI 
CS 
SS 
All 
A12 
A22 
DE 

D19 ( I , J) 


0.0 

1,21 

KM1S1(I , J) / MASS1(I,J) 

KM2S1(I , J) / MASS1(I , J) + E1(I,J) * El(I.J) 
1.0 / ( EA1 (I , J) * FA ) 

BETA1 ( I , J) * PI / 180.0 
DCOS (PITCH) 

DSIN(PITCH) 

CO * CO 
SI * SI 

-EIZ1(I , J) * SS - EIY1(I , J) * CS 
-CO * SI * (EIZ1(I , J) - EIY1(I , J) ) 

EIZ1(I,J) * CS + EIY1 (I , J) * SS 
All * A22 + A12 * A12 
-A 12 * EIY / DE 
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D191(I,J)= A22 * EIY / DE 
D 192(1 ,J)= All * EIY / DE 
D11(I,J) = EIY/GJ1(I,J) 

D12(I , J) = MASS1(I , J)/MASS 
D13(I , J) = D12(I , J) * CO * E1(I , J) / SPAN 
D14(I , J) = D12(I,J) * SI * E1(I , J) / SPAN 
D15(I,J) = D13(I , J) * X * OMEGAN 
D16(I, J) = D14(I , J) * X * OMEGAN 
D17 (I , J) = OMEGAN * D12(I,J) * (KM2S1(I,J) - 
1 KM1S1(I,J)) * (CS - SS) / SPANS 

D18(I,J) = D12(I, J) * (KM1S1(I,J) + KM2S1(I,J)) / SPANS 
X = X + HH 

4251 CONTINUE 

X = SPAN1 / SPAN 

HH = SPAN2 / (100.0* SPAN) 

DO 428 J = 1, 101 

KM1S2( J) = KM1S2 ( J) / MASS2(J) 

KM2S2(J) = KM2S2(J) / MASS2(J) + E2(J) * E2(J) 

EA2(J) = 1.0 / (EA2(J) * FA) 

PITCH = BETA2(J) * PI / 180.0 
CO = DCOS (PITCH) 

SI = DSIN (PITCH) 

CS = CO * CO 
SS = SI * SI 

A12 = -CO * SI * (EIZ2 ( J) - EIY2(J) ) 

All = -EIZ2 ( J) * SS - EIY2(J) * CS 
A22 = EIZ2(J) * CS + EIY2(J) * SS 
DE * All * A22 + A12 * A12 
D29 (J) = -A12 * EIY / DE 

D291(J) = A22 * EIY / DE 
D292(J) = All * EIY / DE 
D21(J) = EIY / GJ2(J) 

D22 ( J) = MASS2 ( J) / MASS 
D23 ( J) = D22 ( J) * CO * E2(J) / SPAN 

D24(J) = D22 ( J) * SI * E2 ( J) / SPAN 

D25 ( J) = D23 (J) * X * OMEGAN 

D26 ( J) = D24(J) * X * OMEGAN 

D27 (J) = OMEGAN * D22(J) * (KM2S2(J) - 

1 KM1S2 ( J) ) * (CS - SS) / SPANS 

D28 (J) = D22(J) * (KM1S2 ( J) + KM2S2(J)) / SPANS 

TEN2(J) = FA * TEN2( J) 

X = X + HH 

428 CONTINUE 

HH1 = SPAN1 / (10.0 * SPAN) 

KH2 = SPAN2 / (50.0 * SPAN) 

STEP = SPAN1 / (20.0 * SPAN) 

H5 = STEP / 24.0 

GO TO (4281, 4282, 4282), NPATH 

4281 BB ( 1) = TEN2 ( 1) 

GO TO 4286 

4282 DO 4284 J = 1, NPATH 

AT ( J) = (9.0 * EA1(J , 1) + 19.0 * EA1(J,2) - 5.0 
1 * EA1(J,3) 4- EA1 ( J ,4) ) * H5 

DO 4283 K = 1,18 

AT (J) = AT ( J) + H5 * ( -EA1(J,K) + 13.0 * EA1(J,K+1) 
+ 13.0 * EAl(J,K+2) - EAl(J,K+4) ) 


1 




4283 CONTINUE 

AT(J) = AT(J) + H5 * ( EA1(J,18) - 5.0 * EA1(J,19) 

1 + 9.0 * EA1(J,20) + 9.0 * EA1(J,21) ) 

4284 CONTINUE 

IF (NPATH .EQ. 3) GO TO 4285 

BB(1) = AT(1) * TEN2(1) / ( AT(1) + AT(2) ) 

BB(2) = AT(2) * TEN2(1) / ( AT(1) + AT(2) ) 

GO TO 4286 

4285 PR = AT(1) * AT(2) + AT(2) * AT(3) + AT(3) * AT(1) 

BB(1) = TEN2(1) * AT (2) * AT(3) / PR 

BB(2) = TEN2(1) * AT(3) * AT(1) / PR 

BB(3) = TEN2(1) * AT(1) * AT(2) /PR 

4286 CONTINUE 
WRITE (22, 380) 

WRITE (22 ,275) ((BB(I) / FA), I = 1, NPATH) 

DO 4290 JJ = 1, NPATH 
X = 0.0 

DO 4287 I = 1,21 

MASS2(I) = MASS1(JJ,I) * H5 * X / MASS 
X = X + STEP 

4287 CONTINUE 
TENl(JJ,21) =0.0 

TEN1(JJ,20) = (9.0 * MASS2 (21) + 19.0 * MASS2(20) 

1 - 5.0 * MASS2(19) + MASS2(18) ) * OMEGAN 

DO 429 I = 2,19 
J = 21 - I 

TEN1(JJ, J) = TEN1(JJ, J+-1) + ( -MASS2(J+2) + 13.0 *(MASS2(J+1) 
1 + MASS2(J)) - MASS2(J-1) ) * OMEGAN 

429 CONTINUE 

TEN1(JJ, 1) = TEN1(JJ,2) + ( MASS2(4) - 5.0 * MASS2(3) 

1 + 19.0 * MASS2(2) + 9.0 * MASS2(1) ) * OMEGAN 

4290 CONTINUE 

IF (ICPL .EQ. 0) ITR = 1 

DO 4291 1=1, NPATH 

DO 4291 J = 1,21 

EIY1(I,J) = TEN1(I,J) + BB(I) 

4291 CONTINUE 
WRITE(22,376) 

DO 4292 I = 1, NPATH 
WRITE (22, 377) I 

WRITE (22, 275) ((EIY1(I, J)/FA) ,J = 1,21) 

4292 CONTINUE 

DO 4318 IK = 1, ITR 
DO 431 J = 1,6 
DO 431 K = 1,6 
XR(J,K) = 0.0 
431 CONTINUE 

CALL TRAMAT(0.0D+00) 

DO 4314 1=1, NPATH 
DO 4310 J = 1,6 
DO 4310 K = 1,6 
XC (J,K) = TFl(I,J,K+6) 

XD(J,K) = TFl(I,J+6,K+6) 

4310 CONTINUE 

CALL SOLUTN (XC , 6 , - 1 , 6 ) 

DO 4311 J = 1,6 
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XC(J,4) = XC(J,4) - FD(I,2) * XC(J,1) 

XC(J, 5) = XC(J,5) - FD ( I , 1 ) * XC(J,1) 

4311 CONTINUE 

CALL MATMUL(6,6 ,6 ,XD,XC ,XE) 

DO 4312 J = 1,6 

DO 4312 K = 1,6 

XC(J,K) = XE(J,K) 

IF (I .EQ. 1) XA(J,K) = XE(J,K) 

IF (I .EQ. 2) XB(J,K) * XE(J,K) 

4312 CONTINUE 

DO 4313 J = 1,6 

XE(1, J) = XE(1, J) - FD ( I , 2 ) * XE(4, J) + FD(I,1) * XE(5,J) 
XE(2, J) = XE(2, J) - FD ( I , 1 ) * XE(6,J) 

XE(3 , J) = XE(3 , J) + FD(I ,2) * XE(6,J) 

4313 CONTINUE 

DO 4314 J = 1,6 

DO 4314 K = 1,6 

XR(J,K) = XR(J,K) + XE(J,K) 

4314 CONTINUE 

CALL SOLUTN (XR , 6 , - 1 , 6 ) 

DO 4315 J = 1,6 

XD( J, 1) = XR(J,6) * TEN2 (1) 

4315 CONTINUE 

DO 4316 1=1, NPATH 

IF (I .EQ. 1) CALL MATMUL(6 ,6 , 1,XA ,XD ,XE) 

IF (I .EQ. 2) CALL MATMUL(6,6,1,XB,XD,XE) 

IF (I .EQ. 3) CALL MATMUL(6 ,6 , 1 ,XC ,XD ,XE) 

BB (I) = XE (6 , 1) 

4316 CONTINUE 

DO 4317 I = 1, NPATH 

DO 4317 J = 1,21 

EIY1CI.J) = TENl(I.J) + BB(I) 

4317 CONTINUE 

WRITE (22, 3 78) IK 
DO 4318 I = 1, NPATH 
WRITE (22, 377) I 

WRITE(22,275) ( (EIY1 (I , J) /FA) , J= 1,21) 

4318 CONTINUE 
WRITE (22, 3 79) 

WRITE (22 ,275) ((TEN2(J)/FA) ,J = 1,101) 

WRITE (22, 205) 

H = H * FACT 
HI = HI * FACT 
H2 = H2 * FACT 
IF (NPATH .LE. 1) GO TO 452 
DO 451 J = 1, NPATH 
DO 451 K = 1,2 
FD(J,K) = FD(J,K)/SPAN 

451 CONTINUE 

452 CONTINUE 

HH1 = SPAN1 / (10.0 * SPAN ) 

HH2 = SPAN 2 / (50.0 * SPAN ) 

IF (I STAGE .EQ. 4) GO TO 471 
IF (I STAGE .NE. 1) GO TO 460 
C 
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noon 


THIS SECTION COMPUTES THE FREQUENCY DETERMINANTS 
IF ISTAGE = 1 


WRITE (22, 200) 

WRITE (22,205) 

WRITE (22, 305) 

455 P = HI * HI 

IF (HI .GT. H2) GO TO S05 
FR = HI/FACT 
F = DET(P) 

WRITE(22,275) FR, HI, F 
HI = H + HI 
GO TO 455 
C 


A-Il 
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THIS SECTION COMPUTES THE NATURAL VIBRATION CHARACTERISTICS 


460 CALL NATFRE(NF) 

IF(IJX .EQ. 0) GO TO 500 
WRITE (22, 200) 

WRITE (22, 205) 

WRITE (22, 3 70) 

WRITE (22, 205) 

DO 461 J a 1, IJK 

FRE = FREQEN(J) / FACT 

HERTZ = FRE / (2.0 * PI) 

WRITE (22, 371) J, FRE, HERTZ 

461 CONTINUE 
WRITE (22, 205) 

IF (I STAGE .EQ. 3) GO TO 505 
ST = SPAN1/ (SPAN * 10.0) 

DO 465 J = 1, 11 

SL1(J) = FLOAT(J-l) * ST 
465 CONTINUE 

ST = SPAN2/ (SPAN * 50.0) 

DO 470 J = 1, 51 

SL2(J) « SL1 ( 11) + FLOAT(J-l) * ST 

470 CONTINUE 
GO TO 473 

471 DO 472 J = 1, NF 
FREQEN(J) = FREQEN(J) * FACT 

472 CONTINUE 
IJK = NF 

473 DO 495 IJ = l.IJK 

J1 = IJ 

PP = FREQEN(IJ) 

P = PP * PP 

CALL SHAPES (P ,W1 ,W2 , VI , V2 ,PHI 1 , PHI2 ,U ,UC , PSIC , ANUC) 

AMAX = Wl(l,l) 

DO 475 I = l.NPATH 
DO 475 J = 1, 11 

IF(DABS(AMAX) .LT. DABS (W1(I , J) ) ) AMAX = W1(I,J) 

IF (DABS (AMAX) .LT. DABS (VI ( I , J) ) ) AMAX = V1(I,J) 

IF (DABS (AMAX) .LT. DABS(PHI1(I , J) ) ) AMAX = PHI1(I,J) 
475 CONTINUE 

DO 480 J = 1, 51 

IF (DABS (AMAX) .LT. DABS(W2(J))) AMAX = W2(J) 

IF (DABS (AMAX) .LT. DABS(V2(J))) AMAX = V2(J) 

IF (DABS (AMAX) .LT. DABS (PHI2 ( J) ) ) AMAX = PHI2(J) 

480 CONTINUE 

DO 485 I = 1, NPATH 

DO 485 J = 1, 11 
U(I) = U ( I ) /AMAX 
Wl(I.J) = W 1 ( I , J ) /AMAX 
V1(I,J) = VI (I ,J)/AMAX 
PHU (I , J) = PHI1(I , J)/AMAX 
485 CONTINUE 

DO 490 J = 1,51 

W2(J) = W2 ( J) /AMAX 


A-12 



490 


496 


495 

500 

505 

C 


V2(J) = V2(J)/AMAX 
PHI2(J) * PHI2(J)/AMAX 
CONTINUE 

FRE = PP / FACT 
HERTZ = FRE / (2.0 * PI) 

CPM = HERTZ * 60.0 
CALL PLOT ( 1 , W2, Wl) 

CALL PLOT (2, V2, VI) 

CALL PLOT (3, PHI2, PHI1) 

DO 496 J * 1, NPATH 
U(J) = U(J) / AMAX 
WRITE (22, 372) J, U(J) 

CONTINUE 

PSIC = PSIC / AMAX 
WRITE(22 ,373) PSIC 
ANUC = ANUC / AMAX 
WRITE (22, 374) ANUC 
CONTINUE 

IF(IJK .LT. NF) WRITE (22, 320) IJK 
IF (I STAGE .EQ. 1) WRITE(22,205) 
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FORMATS 


100 FORMAT (5 15) 

105 FORMAT (5E 14. 7) 
110 FORMAT (3A1) 

200 FORMAT (1H1) 

205 FORMAT (//2X,' 


210 FORMAT (//5X, 'NUMBER OF FREQUENCIES REQUIRED = ’,15) 

215 FORMAT (//5X, 'FREQUENCY INCREMENTS (RAD/SEC) = ’,E14.7) 

220 F0RMAT(//5X, 'ENDING FREQUENCY (RAD/SEC) = ’,E14.7) 

225 FORMAT (//5X, 'RADIUS OF THE ROTOR (INCHES) = ’,E14.7) 

230 FORMAT (//5X, 'LENGTH OF INBOARD SEGMENTS (IN)= ’.E14.7) 

235 F0RMAT(//5X, 'LENGTH OF THE BLADE (INCHES) = ’.E14.7) 

236 FORMAT (//5X, 'SEMI -CHORD (INCHES) = ',E14.7) 

240 F0RMAT(//5X, 'ROTATIONAL SPEED (RPM) = ’.E14.7) 

245 FORMAT (//5X, 'NUMBER OF DATA POINTS FOR INBOARD SEGMENTS = ',15) 
250 FORMAT (//5X, 'NUMBER OF DATA POINTS FOR BLADE = ',15) 

255 FORMAT (//5X, 'PROPERTIES OF THE FIRST LOAD PATH') 


260 FORMAT (//5X, 'PROPERTIES OF THE SECOND LOAD PATH') 

265 F0RMAT(//5X, 'PROPERTIES OF THE THIRD LOAD PATH') 

270 F0RMAT(//5X, 'DATA POINT LOCATIONS IN INCHES’) 

275 FORMAT (4(6X,El4.7)) 

280 FORMAT (//5X, 'MASS PER UNIT LENGTH (LB-SEC .**2/IN**2) ' ) 

285 F0RMAT(//5X, 'FLAPWISE BENDING STIFFNESS (LB-IN**2)’) 

286 FORMAT (//5X, 'AXIAL STIFFNESS (LB) ’) 

290 F0RMAT(//5X, 'PROPERTIES OF THE BLADE’) 

295 F0RMAT(//5X,' INTERPOLATED VALUES FOR INBOARD SEGMENTS, 21 ’, 

1 'EQUIDISTANT VALUES') 

300 FORMAT (//5X, 'INTERPOLATED VALUES FOR THE BLADE, 101 ’, 

1 'EQUIDISTANT VALUES') 

305 F0RMAT(//1X, 'FREQUENCY (RAD/SEC) NON-DIMENSIONAL FREQUENCY ', 
1' DETERMINANT') 

310 F0RMAT(//5X, 'STARTING FREQUENCY (RAD/SEC) = ’,E14.7) 

315 FORMAT (//5X, 'NUMBER OF LOAD PATHS = ',15) 

320 F0RMAT(//5X, 'NUMBER OF FREQUENCIES DETECTED WITHIN RANGE', 

1/5X, ' = ’,15) 

321 F0RMAT(//5X, 'Y, Z DISTANCES BETWEEN LOAD PATH NO: ',12, 

1/5X, ' AND THE BLADE (IN) ARE = ' ,F6 .2 ,5X,F6 .2) 

343 F0RMAT(//2X, 'AXIAL STIFFNESS OF THE BLADE (LB) IS ') 

347 FORMAT (//5X, 'EXECUTION STAGE REQUIRED IS = ',15) 

348 FORMAT (//5X, 'PLOT INDEX IS = *,I5) 

350 FORMAT (//5X, 'CHORDWISE BENDING STIFFNESS (LB- IN**2) ' ) 

352 FORMAT (//5X, 'TORSIONAL STIFFNESS (LB -IN**2) ’ ) 

354 FORMAT (//5X, 'DISTANCE BETWEEN MASS AND ELASTIC AXIS(IN) =') 

356 F0RMAT(//5X, 'TWIST OF THE BLADE INCLUDING THE COLLECTIVE’, 

1/5X, ' (DEGREES) = ’) 

358 F0RMAT(//5X, 'MASS MOMENT OF INERTIA ABOUT THE CHORD ( LB - SEC**2) ' ) 
360 FORMAT (//5X, 'MASS MOMENT OF INERTIA ABOUT AN AXIS ’, 

1/5X, 'PERPENDICULAR TO THE CHORD THROUGH THE CENTER OF ', 

1 /5X , ' GRAVITY (LB-SEC**2 ) ' ) 

370 FORMAT (/5X, 'THE NATURAL FREQUENCIES OF THE SYSTEM ARE:', 

1/17X, ' RAD / SEC: ', ' HERTZ: ’) 

371 F0RMAT(5X, 15, 5X, F11.4, 5X, F11.4) 



372 FORMAT (//5X, 'AXIAL DISPLACEMENT OF LOAD PATH # ' ,I3,2X, 
1/SX, 'AT THE CLEVIS = ’.E11.4) 

373 FORMAT(//5X, 'FLAPWISE BENDING SLOPE AT THE CLEVIS = \E11.4) 

374 FORMAT (//5X, 'CHORDWISE BENDING SLOPE AT THE CLEVIS = ’,E11.4) 

375 FORMAT(//5X, 'NUMBER OF TENSION ITERATIONS = ’,15) 

376 FORMAT C//5X, 'THE FOLLOWING ARE THE STARTING TENSIONS') 

377 FORMAT (//5X, 'THE TENSION COEFFICIENTS IN LOAD' 

1/5X, 'PATH #' ,15,' ARE') 

378 FORMAT(//5X, 'THE TENSIONS AFTER ITERATION #’,I5, ' ARE’) 

379 FORMAT(//5X, 'THE TENSIONS IN THE BLADE ARE’) 

380 FORMAT (//5X, 'STARTING TENSIONS IN THE LOAD PATHS AT THE' 

1 /5X, 'CLEVIS ARE') 

STOP 

END 

C 
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FUNCTION DET(P) 

IMPLICIT REAL*8 (A - H, 0 - Z) 

DIMENSION TF1(3, 12, 12), TF2(12,12), A(6,6), B(6,6), 

1 C(6,6) , D(6,6) , E(6 ,6) , R(6,6), FD(3,2) 

C0MM0N/X4/NPATH , ICPL 
COMMON/X9/FD 
COMMON/XIO/TFl, TF2 
COMMON /XS 2 /DETER 
CALL TRAMAT(P) 

DO 10 I = 1, 6 
DO 10 J = 1, 6 

A(I , J) = TF2(I+6, J) 

B(I , J) = TF2(I+6, J+6) 

R(I, J) « 0.0 
10 CONTINUE 

DO 20 1 = 1, NPATH 

DO 15 J = 1, 6 
DO 15 K = 1, 6 
C ( J,K) = TF1(I , J, K+6) 

D(J,K) = TF1(I , J+6 , K+6) 

15 CONTINUE 

CALL SOLUTN (C , 6, -1, 6) 

DO 16 J = 1, 6 

C(J,4) = C (J ,4) - FD ( I , 2 ) * C(J,1) 

C(J,5) a C(J,5) - FD(I,1) * C(J,1) 

16 CONTINUE 

17 CALL MATMUL(6 , 6, 6, D, C, E) 

DO 18 J = 1,6 

E(l, J) = E(1 , J) - FD(I ,2) * E(4, J) + FD(I,1) * E(5,J) 
E(2,J) = E(2 , J) - FD ( I , 1 ) * E(6,J) 

E(3 , J) = E(3 , J) + FD ( I , 2 ) * E(6,J) 

18 CONTINUE 

DO 20 J = 1, 6 
DO 20 K a 1, 6 
R(J,K) = R(J,K) + E(J,K) 

20 CONTINUE 

CALL MATMUL(6, 6, 6, B, R, D) 

DO 25 J = l, 6 
DO 25 K = 1, 6 

B(J,K) a D(J,K) + A(J,K) 

25 CONTINUE 

CALL SOLUTN (B , 6, -1, 6) 

DET = DETER 

RETURN 

END 
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SUBROUTINE INTPOL(N, A, H) 

C 

C 

C THIS SUBROUTINE INTERPOLATES FOR THE REQUIRED VALUES 

C - - 

C 

IMPLICIT REAL*8 (A • H, 0 • Z) 

DIMENSION A (101), STA(lOl), TABLE(101,1) , B(101) 

C0MM0N/X3/STA , NS 
NN = N - 1 
A(N) = A(NS) 

NM1 = NS -1 
DO 20 I = 1, NM1 

20 TABLE (I , 1) = (A(I+1) - A(I) )/(STA(I+l) - STA(I)) 

XARG = H 
DO 35 I = 2, NN 
DO 25 J = 1, NS 

IF (J .EQ. NS .OR. XARG .LE. STA(J)) GO TO 30 
25 CONTINUE 

30 MAX = J 

IF (MAX .LE. 2) MAX = 2 

ISUB a MAX - 1 

YEST = TABLE ( ISUB, 1) 

B(I) = YEST * (XARG - STA(ISUB)) + A(ISUB) 

35 XARG = XARG + H 
DO 40 J = 2, NN 
40 A(J) = B (J) 

RETURN 

END 


A-17 
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SUBROUTINE MATMUL(L, M, N, A, B, C) 
C 

C 

C MATRIX MULTIPLICATION 

C 

C 

IMPLICIT REAL*8 (A - H, 0 - Z) 

DIMENSION A(6 ,6) , B(6,6), C(6,6) 

DO 10 I = 1, L 
DO 10 J = 1, N 
C(I,J) = 0.0 
DO 10 K = 1, M 

C(I, J) = C(I, J) + A(I,K) * B(K,J) 

10 CONTINUE 
RETURN 
END 
C 
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SUBROUTINE MROOT(SDT, H, FRE, IJ, ICOUNT) 


THIS SUBROUTINE CALCULATES THE MISSING FREQUENCIES, IF ANY 


10 

15 


20 


30 


40 

50 


IMPLICIT REAL*8 (A - H, 0 - Z) 

DIMENSION SDT(2,2000) , FREMIS(IO), FRE(IO) 

C0MM0N/X2/FACT 

INDEX = 0 

N = ICOUNT - 2 

IJ = 0 

DO 10 I = 1, N 
A = SDT(2,I) 

B = SDT (2,1+1) 

C = SDT (2, 1+2) 

D = A * B 

IF(D .LT. 0.0) GO TO 10 
D = B * C 

IF(D .LT. 0.0) GO TO 10 
D = DABS (A) - DABS(B) 

IF(D .LT. 0.0) GO TO 10 
D = DABS(C) - DABS(B) 

IF(D .LT. 0.0) GO TO 10 
INDEX = INDEX + 1 
IF ( INDEX .GT. 5) GO TO 15 
FREMISC INDEX) * SDT (1,1) 

CONTINUE 

IF ( INDEX .EQ. 0) GO TO 90 
CONTINUE 

DO 60 I = 1, INDEX 
NS = 1 

HH = H/ (FLOAT (NS) * 10.0) 

PP = FREMIS( INDEX) 

P = PP * PP 
F = DET(P) 

F = DSIGN(1.0D+00,F) 

A = FREMIS( INDEX) + 2.0 * H 
PP = PP + HH 

IF (PP .GT. A) GO TO 50 
P = PP * PP 
G = DET(P) 

G = DSIGN( 1 . OD+OO ,G) 

IF (F*G .GT. 0.0) GO TO 40 
IJ = IJ + 1 
P = (PP-HH) * (PP-HH) 

C = DET(P) 

P = PP * PP 
D = DET(P) 

FRE (IJ) = PP - D * HH/(D - C) 
IF((IJ/2*2-IJ) .EQ. 0) GO TO 60 
F = G 

GO TO 30 
NS = NS * 10 
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IF (NS .GT. 100) GO TO 70 
GO TO 20 

60 CONTINUE 
70 WHITE (22, 100) 

WRITE (22, 110) 

DO 80 I * 1, INDEX 

A = FREMIS(I)/FACT 
B = A + 2.0 * H/FACT 
WRITE (22, 120) A, B 
80 CONTINUE 
90 CONTINUE 
100 FORMAT (1H1) 

110 FORMAT (10X,’ CHECK FOR TWO FREQUENCIES BETWEEN EACH OF THE 
1' ,/10X, 'FOLLOWING SETS. OTHERWISE, THEY ARE MISSING') 

120 FORMAT (2 ( 10X,E14 . 7) ) 

RETURN 

END 

C 
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SUBROUTINE NATFRE(N) 

c 

C 

C THIS SUBROUTINE SCANS THE FREQUENCY DETERMINANT WITH RESPECT TO 
C THE FREQUENCY UNTIL THE SPECIFIED NUMBER OF SIGN CHANGES ARE 
C DETECTED, STARTING FROM ZERO FREQUENCY. USES THE "DET" ROUTINE 

C 

C 

IMPLICIT REAL*8 (A - H, 0 ■ Z) 

DIMENSION FREQEN(IO), JKL(IO), SDT(2,2000), FRE(IO), STO(20) 
COMMON/X1/FREQEN , H, HI, H2, IJK 
IJK =0 
IJ =0 
ICOUNT = 1 
DO 5 J = 1, N 
5 JKL(J) = 0 
PP = HI 
P = PP * PP 
F = DET(P) 

SDTC1, ICOUNT) = PP 
SDT (2, ICOUNT) = F 
IF (DABS (F) .GT. O.OOOI) GO TO 10 
IJK = IJK + 1 

JKL(IJK) = 1 

FREQEN(IJK) = PP 

PP * PP + H 

p s PP * pp 

F = DET(P) 

ICOUNT = ICOUNT + 1 
SDT (1, ICOUNT) = PP 
SDT (2, ICOUNT) = F 
10 F = DSIGN ( 1 . OD+OO ,F) 

15 PP = PP + H 

IF (PP .GT. H2) GO TO 55 
P = PP * PP 
G = DET(P) 

ICOUNT = ICOUNT + 1 

SDT (1, ICOUNT) = PP 

SDT(2, ICOUNT) = G 

IF (DABS (G) .GT. 0.0001) GO TO 20 

IJK = IJK + 1 

JKL(IJK) = 1 

FREQEN(IJK) = PP 

IFUJK .EQ. N) GO TO 55 

PP = PP + H 

P = PP * PP 

F = DET(P) 

ICOUNT = ICOUNT + 1 
SDT (1, ICOUNT) = PP 
SDT (2, ICOUNT) = F 
F = DSIGN(1. OD+OO, F) 

GO TO 15 

20 G = DSIGN(1. OD+OO, G) 

IF (F*G .GT. 0.0) GO TO 25 


IJK = IJK + 1 
FREQEN(IJK) = PP - H 
IF (IJK .EQ. N) GO TO 55 
25 F = G 

GO TO 15 

55 IF (IJK .EQ. 0) GO TO 65 
HS = H/10.0 
DO 50 J = 1, IJK 

IF(JKL(J) .EQ. 1) GO TO 50 
PP = FREQEN(J) 

P = PP * PP 
F = DET(P) 

F = DSIGN(1. OD+OO.F) 

35 PP = PP + HS 

P a PP * PP 
G = DET(P) 

IF (DABS (G) .GT. 0.0001) GO TO 40 
JKL(J) = 1 
FREQEN(J) = PP 
GO TO 50 

40 G = DS IGN(1.0D+00,G) 

IF(F*G .GT. 0.0) GO TO 45 
FREQEN(J) = PP - HS 
GO TO 50 
45 F = G 

GO TO 35 
50 CONTINUE 

DO 60 J = 1, IJK 

IF(JKL(J) .EQ. 1) GO TO 60 
PP = FREQEN(J) 

P = PP * PP 
F = DET(P) 

PP = PP + HS 
P = PP * PP 
G = DET(P) 

DIFF = G - F 

FREQEN(J) = PP - G*HS/DIFF 
60 CONTINUE 
65 CONTINUE 

CALL MROOT ( SDT , H, FRE, IJ, ICOUNT) 
IF(IJ .EQ. 0) GO TO 100 
N2 = IJK + IJ 
IF (IJK .EQ. 0) GO TO 75 
DO 70 I = 1, IJK 
70 STO(I) = FREQEN(I) 

75 N1 = IJK + 1 

DO 80 I = Nl, N2 
80 STO(I) = FREQEN(I-IJK) 

DO 90 1=2, N2 
DO 90 J = I, N2 

IF(ST0(I-1) - STO(J) ) 90, 90, 85 
85 STORE = STO(J) 

STO(J) = STO(I-l) 

STO(I-l) = STORE 
90 CONTINUE 
IJK = N2 



IF(N2 .GT. N) IJK = N 
DO 95 I = 1, IJK 
95 FREQEN(I) = STO(I) 

100 CONTINUE 
RETURN 
END 
C 
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SUBROUTINE PLOT(NN, C, B) 


PRINTS AND PLOTS NATURAL VIBRATION CHARACTERISTICS 


50 


60 


80 


81 


120 

130 

140 


IMPLICIT REAL*8 (A - H, 0 - Z) 

REAL*8 LINE (61) 

DIMENSION C(51), B(3,ll), SLl(ll), SL2(51), A(61), SL(61) 
COMMON /X4/NPATH , ICPL 

C0MM0N/X6 / SL 1 , SL2 , CPM , FRE , HERTZ , BLANK .DOT, STAR , J 1 , I PLOT 
WRITE(22,10) 

WRITE (22, 20) 

WRITE (22, 30) Jl, FRE, HERTZ, CPM 
WRITE (22 ,20) 

IF(NN .EQ. 1)WRITE(22,39) 

IF(NN .EQ. 2)WRITE(22,40) 

IF(NN .EQ. 3)WRITE(22,41) 

DO 50 J = 1, 11 
SL( J) = SL1(J) 

A( J) = B ( 1 , J) 

CONTINUE 

DO 60 J = 2, 51 

SL( J+10) = SL2 ( J) 

A( J+10) = C(J) 

CONTINUE 
WRITE (22, 70) 

DO 80 J = 1, 11 

WRITE (22 ,90) SL(J) ,A(J) ,SL(J+21) ,A(J+21) ,SL(J+42) ,A(J+42) 

CONTINUE 

DO 81 J = 12,19 

WRITE(22,90) SL(J) ,A(J) ,SL(J+21) ,A(J+21) ,SL(J+42) ,A(J+42) 
CONTINUE 

WRITE (22 ,90) SL(20) , A(20) , SL(41), A(41) 

WRITE(22,90) SL(21) , A(21), SL(42), A(42) 

WRITE (22, 20) 

IF(NPATH .LE. 1) GO TO 140 
WRITE (22, 10) 

WRITE (22, 20) 

DO 130 I = 2, NPATH 
IF(NN .EQ. 1) WRITE (22 , 100) I 
IF(NN .EQ. 2) WRITE(22 , 101) I 
IF(NN .EQ. 3) WRITE (22 , 102) I 
WRITE (22, 70) 

DO 120 J = 1, 3 

WRITE (22 ,90) SL1(J) ,B(I,J) ,SLl(J+4) ,B(I,J+4) ,SLl(J+8) , 

1 B (I , J+8) 

CONTINUE 

WRITE(22, 90) SL1(4), B(I,4), SL1(8), B(I,8) 

CONTINUE 
WRITE (22, 20) 

WRITE (22, 10) 

IF(IPLOT .EQ. 0) GO TO 180 
DO 150 J = 1, 61 


LINE(J) = DOT 
150 CONTINUE 

J = 30.0 * (A(l) +1.0) + 1.5 
LINE(J) = STAR 

WRITE(22,110) (LINE(J) , J = 1,61) 

DO 160 J s i, 61 
LINE ( J) = BLANK 
160 CONTINUE 

LINE (31) = DOT 
DO 170 JJ = 2, 61 

J = 30.0 * (A(JJ) + 1.0) + 1.5 
IF(J .GT. 61) GO TO 170 
LINE (J) = STAR 

WRITE(22, 190) (LINE(JV), JV = 1,61) 

LINE ( J) = BLANK 
LINE (31) = DOT 
170 CONTINUE 
180 CONTINUE 
10 FORMAT (1H1) 

20 FORMAT (/2X,’ 

30 FORMAT (//5X,' MODE NUMBER = ' ,I2,8X, ’FREQ. RAD / SEC = ’,F10.4,8X, 

1/5X, 'FREQ. HERTZ = ' ,F10 .4,8X, ’FREQ. CYCLES/MINUTE = ’.F10.4) 

39 FORMAT (//15X,* FLAPWISE DEFLECTION/SEMI-CHORD ') 

40 FORMAT (//15X,' CHORDWISE DEFLECTION/SEMI -CHORD ') 

41 F0RMAT(//15X, ' TORSIONAL DEFLECTION (RAD) ’ ) 

70 FORMAT (//4X,’STA X/L' ,3X, ’DEFLN' ,9X, 'STA X/L' ,3X, ’DEFLN’ ,9X, 
l'STA X/L' ,3X, 'DEFLN' ) 

90 F0RMAT(//3(2X,F8.4,3X,E11.4)) 

100 FORMAT (//5X, ’FLAPWISE DEFLECTION/SEMI -CHORD LOAD PATH # ’,15) 

101 F0RMAT(//5X, ’CHORDWISE DEFLECTION/SEMI -CHORD LOAD PATH # ’ , 15 ) 

102 FORMAT (//5X, 'TORSIONAL DEFLECTION (RAD) LOAD PATH # ’,15) 

110 FORMAT (//4X, 6 1A1) 

190 FORMAT (4X, 6 1A1) 

1101 RETURN 
END 
C 


/ 
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FUNCTION RUNGE ( Y , F, J, M, HH) 

Mr k4r ^tt r ^ ' \ ^ J r &r & >4r&r& ' J: - ir^J r Jr^* - J r k - Jrir{rX - JrX - i r i; - A ■■ !. ■ ' A - A - » - A - A ~ A - A - - A - A ■ A - A • k&M rMc lriririr it 


FOURTH ORDER RUNGE -KUTTA METHOD 


IMPLICIT REAL*8 (A - H, 0 - Z) 

INTEGER RUNGE 

DIMENSION Y(12), F(12), PHI(12), SAVEY(12) 

NN = 12 
M ■ M + 1 

GO TO (5 , 10, 20, 30, 40), M 
5 RUNGE = 1 
RETURN 

10 DO 15 JJ = 1, NN 

SAVEY(JJ) = Y(JJ) 

PHI ( JJ) = F( JJ) 

Y(JJ) = SAVEY(JJ) + HH * F(JJ) / 2.0 
15 CONTINUE 
J = J + 1 
RUNGE = 1 
RETURN 

20 DO 25 JJ = 1, NN 

PHI(JJ) = PHI(JJ) + 2.0 * F ( JJ) 

Y ( JJ) = SAVEY(JJ) + HH * F(JJ) / 2.0 
25 CONTINUE 
RUNGE = 1 
RETURN 

30 DO 35 JJ = 1, NN 

PHI ( JJ) = PHI ( JJ) + 2.0 * F(JJ) 

Y ( JJ) = SAVEY(JJ) + HH * F(JJ) 

35 CONTINUE 
J=J+1 
RUNGE = 1 
RETURN 

40 DO 45 JJ = 1, NN 

45 Y ( JJ) = SAVEY(JJ) + (PHI(JJ) + F(JJ))* HH / 6.0 

M = 0 
RUNGE = 0 
RETURN 
END 
C 
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SUBROUTINE SHAPES(P,W1,W2,V1,V2,PHI1,PHI2,U,UC,PSIC,ANUC) 

C j W ^rthfr ^ T H r A^ fr ^^^ - A^A - tVA ' A - A - A ^ A - ^ - A^^^ - A n ^ - A - A - A - ^^ 

C 

c - - 

C THIS SUBROUTINE COMPUTES THE NATURAL MODE SHAPES 

C 

C 

IMPLICIT REAL*8 (A - H, 0 - Z) 

DIMENSION Wl(3,ll) ,W2(51) ,V1(3,11) , V2(51), 

1 PHI1(3,11), PHI2(S1) ,U(3) ,TF1(3 , 12 , 12) , 

1TF2(12,12) ,TT1(3,11,3,12) ,TT2 (51 ,3 , 12) , 
1A(6,6),B(6,6),C(6,6),D(6,6),E(6,6), 

IFD(3,2) ,R(6,6),BB(12) ,TEM(12,12) , G(5,5) 

C0MM0N/X4/NPATH , ICPL 
COMMON /X8 /BB , IND 
COMMON/X9/FD 
COMMON /X 10 /TF 1 , TF2 
COMMON/X12/CON 
COMMON/XS3/TT1 , TT2 
CALL TRAMAT(P) 

DO 5 I = 1,12 
DO 5 J = 1,12 
TEM(I , J) = TF2(I,J) 

5 CONTINUE 

DO 10 I = 1,6 

DO 10 J = 1,6 

A(I, J) = TF2(I+6,J) 

B(I, J) ■» TF2CI+6 , J+6) 

R(I , J) * 0.0 
10 CONTINUE 

DO 20 I = 1, NPATH 
DO 15 J = 1,6 

DO 15 K = 1,6 

C ( J ,K) = TF1(I , J ,K+6) 

D(J,K) = TF1 (I , J+6 ,K+6) 

15 CONTINUE 

CALL SOLUTN (C , 6 , - 1 , 6 ) 

DO 16 J = 1,6 

C(J,4) * C(J,4) - FD ( 1 , 2 ) * C( J , 1) 

C(J,5) = C(J,5) - FD ( 1 , 1 ) * C ( J , 1) 

16 CONTINUE 

17 CALL MATMUL(6 ,6,6,D,C,E) 

DO 18 J = 1,6 

E(l, J) a E ( 1 , J) - FD ( I , 2 ) * E (4 , J)+ FD(I,1) * E(5,J) 

E(2 , J) = E (2 , J) - FD ( I , 1 ) * E(6 , J) 

E (3 , J) = E (3 , J) + FD(I ,2) * E(6,J) 

18 CONTINUE 

DO 20 J = 1,6 
DO 20 K = 1,6 
R(J,K) a R(J,K) + E (J ,K) 

20 CONTINUE 

CALL MATMUL (6,6,6,B,R,D) 

DO 25 J = 1,6 
DO 25 K = 1,6 
E (J,K) = D(J,K) + A(J,K) 
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25 CONTINUE 

CALL SOLUTN (TEM ,12,-1,12) 

DO 30 J = 1,6 

DO 30 K = 1,6 

D(J,K) = TEM(J,K) 

30 CONTINUE 

CALL MATMUL(6,6,6,E,D,A) 

WRITE(24,111) ( (A(I , J) , J = 1,6), I = 1,6) 

111 FORMAT (/ (6 (Ell. 4, IX) )) 

IFflCPL .NE. 0) GO TO 310 

C(l,l) = 0.0 

C(2,l) = 1.0 

DO 32 K = 1,4 

DO 31 J = 1,4 

G(J,K) = A(J,K+2) 

31 CONTINUE 
G(K,5) = -A(K,2) 

32 CONTINUE 

CALL SOLUTN (G , 4 , 1 , 5 ) 

DO 33 J = 3,6 
C(J,1) = BB ( J-2) 

33 CONTINUE 
GO TO 375 

310 DO 36 K = 1,5 
KK = K 

IF(K .GE. 2) KK = K + 1 
DO 35 J = 1,5 
D(J,K) = A(J,KK) 

35 CONTINUE 

D(K,6) = -A(K,2) 

36 CONTINUE 

CALL SOLUTN (D,5,l,6) 

C(l,l) = BB ( 1) 

C (2 , 1) = 1.0 
DO 37 J = 3,6 
C(J,1) = BB(J-l) 

37 CONTINUE 

375 CONTINUE 

DO 376 I = 1,6 
A(1 , 1) = 0.0 
B(I,1) = 0.0 
DO 376 K = 1,6 

A(I , 1) = A(I , 1) + TEM(I,K) * C (K, 1) 

B (I , 1) = B(I , 1) + TEM (1+6 ,K) * C(K,1) 

376 CONTINUE 

DO 39 I = 1,51 
DO 39 J = 1,3 
STO = 0.0 
DO 38 K = 1,6 

STO = STO + TT2 (I , J ,K) * A(K,1) + TT2(I,J,K+6) * B(K,1) 

38 CONTINUE 
PHI2(I) = STO 

IF ( J .EQ. 1) W2(I) = STO * CON 
IF ( J .EQ. 2) V2(I) = STO * CON 

39 CONTINUE 

DO 48 J = 2,6 
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C(J,1) = A(J , 1) 

48 CONTINUE 

DO 60 I = l.NPATH 

C(1,1)=A(1,1)-FD(I,2)*A(4,1)-FD(I,1)*A(5,1) 

DO 50 J = 1,6 
DO 50 K = 1,6 
B(J,K) = TFi(I,J,K+6) 

50 CONTINUE 

CALL SOLUTN (B , 6 , - 1 , 6 ) 

CALL MATMUL(6 ,6 , 1,B ,C ,D) 

DO 55 J = 1,11 
W1(I,J) = 0.0 
V1(I,J) = 0.0 
PHI1(I,J) = 0.0 
DO 55 K = 1,6 

W1(I,J) = W1(I,J) + ( TTl(I,J,l,K+6) * D(K, 1) ) * CON 
V1(I,J) = V1(I,J) + ( TT 1 ( I , J , 2 , K+6 ) * D(K, 1) ) * CON 
PHI1(I , J) = PHIKI.J) + TTl(I,J,3,K+6) * D(K,1) 

55 CONTINUE 

U(I) = C(l, 1) 

60 CONTINUE 

UC = A(l, 1) 

PSIC = A(4 , 1) 

ANUC = A(5 , 1) 

RETURN 

END 

C 
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c 

c 

c 


50 

55 


60 


65 


70 


75 

80 


85 


SUBROUTINE SOLUTN(A, N, INDIC, NRC) 


IMPLICIT REAL*8 (A - H, 0 • Z) 

DIMENSION A(NRC,NRC) ,X(12) ,IR0W(12) ,JC0L(12) ,J0RD(12) ,Y(12) 

C0MM0N/XS2/DETER 

C0MM0N/X8/X, IND 

IND = 0 

MAX = N 

IF (INDIC .GE. 0) MAX = N + 1 
DETER =1.0 
DO 80 K = 1, N 
KM1 = K - 1 
PIVOT =0.0 
DO 60 I = 1, N 
DO 60 J = 1, N 

IF (K .EQ. 1) GO TO 55 
DO 50 ISCAN = 1, KM1 
DO 50 JSCAN = 1, KM1 

IF (I .EQ. IROW(ISCAN) ) GO TO 60 
IF(J .EQ. JCOL( JSCAN) ) GO TO 60 
CONTINUE 

IF (DABS (A ( I, J) ) .LE. DABS(PIVOT)) GO TO 60 
PIVOT = A(I , J) 

IROW(K) = I 
JCOL(K) = J 
CONTINUE 

IF (DABS (PIVOT) .GT. 0.1E-20) GO TO 65 
DETER =0.0 
IND = 1 
RETURN 

IROWK = IROW(K) 

JCOLK = JCOL(K) 

DETER = DETER * PIVOT 

DO 70 J = 1, MAX 

A(IROWK , J) = A ( IROWK , J ) / P I VOT 

A (IROWK, JCOLK) = l.O/PIVOT 

DO 80 I = 1, N 

AIJCK = A(I, JCOLK) 

IF (I .EQ. IROWK) GO TO 80 
A (I, JCOLK) = -AIJCK/PIVOT 
DO 75 J = 1, MAX 

IF ( J .NE. JCOLK) A(I,J) = A(I,J) - AIJCK * A(IROWK,J) 
CONTINUE 
DO 85 I = 1, N 
IROWI = IROW(I) 

JCOLI = JCOL(I) 

JORD( IROWI) = JCOLI 

IF (INDIC .GE. 0) X(JCOLI) = A(IROWI,MAX) 

CONTINUE 
INTCH = 0 
NM1 = N-l 
DO 90 I = 1, NM1 
IP1 =1+1 
DO 90 J = IP1, N 
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IF(JORD(J) .GE. JORD(I) ) GO TO 90 
JTEMP = JORD(J) 

JORD(J) = JORD(I) 

JORD(I) = JTEMP 
INTCH = INTCH + 1 
90 CONTINUE 

IF (INTCH/ 2*2 .NE. INTCH) DETER = -DETER 

IFCINDIC .LE. 0) GO TO 94 

RETURN 

94 DO 100 J » 1, N 

DO 95 I = 1, N 
IROWI = IROW(I) 

JCOLI = JCOL(I) 

Y(JCOLI) = A( IROWI , J) 

95 CONTINUE 

DO 100 I = 1, N 
A(I,J) = Y(I) 

100 CONTINUE 

DO 110 I = 1, N 
DO 105 J = 1, N 
IROWJ = IROW(J) 

JCOLJ = JCOL(J) 

Y( IROWJ) = A (I, JCOLJ) 

105 CONTINUE 

DO 110 J = 1, N 
A(I , J) = Y (J) 

110 CONTINUE 

RETURN 
END 
C 
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SUBROUTINE TRAMAT(P) 

C 

C - 

C THIS SUBROUTINE COMPUTES THE TRANSFER MATRIX 

C 

C 

IMPLICIT REAL*8 (A - H, 0 - Z) 

INTEGER RUNGE 

DIMENSION TF1 (3,12,12) , TF2(12,12), TT1(3 , 11,3 , 12) , 
1TT2(51,3 , 12) ,T(12) ,V(12) , TEN1(3,21), TEN2(101), EA1(3,21), 
1EA2(101) ,D11(3 ,21) ,D21(101) ,D12(3,21) ,D22(101) ,D13(3,21) , 
1D23(101) ,D14(3,21) ,D24(101) ,D1S(3,21) ,D25(101) ,D16(3,21) , 
1D26 (101),D17(3,21),D27(101),D18(3,21) ,D28 (101),D19(3,21), 
1029(101) ,D191(3 ,21) ,D291(101) ,D192(3,21) ,D292(101) 
C0MM0N/X4/NPATH , ICPL 

C0MM0N/X5 /TEN 1 , TEN2 ,EA1 ,EA2 ,D11 ,D21 ,D12 ,D22 ,D13 , 

1D23 ,D14 ,D24 ,D15 ,D25 ,D16 ,D26 ,D17 ,D27 ,D18 ,D28 ,D19 , 

1D29 ,D191 ,D291 ,D192 ,D292 .OMEGAN 
C0MM0N/X7 /HH1 ,HH2 
COMMON/XIO/TFl, TF2 
C0MM0N/XS3/TT1 , TT2 
HH = HH1 
PP = P + OMEGAN 
DO 35 I = 1, NPATH 
DO 35 J = 1, 12 
DO 10 K = 1, 12 
V(K) = 0.0 
10 CONTINUE 

V(J) = 1.0 
M = 0 
NJ = 1 

DO 25 L = 1, 10 

15 K = RUNGE (V, T, NJ, M, HH) 

IF (K .NE. 1) GO TO 20 
T( 1) = EA1(I,NJ)*V(12) 

T(2) = V(4) 

T(3) = V(5 ) 

T(4) = D19 (I ,NJ) * V(8) + D191(I,NJ) * V(9) 

T(5) = D192 ( I ,NJ) * V(8) - D19(1,NJ) * V(9) 

T(6) = D11(I,NJ) * V(7 ) 

T(7) = -P * D13(I,NJ) * V(2) + PP * D14(I ,NJ) * V(3) 

1 + D15 (I ,NJ) * V(4) - D16 (I ,NJ) * V(5) 

1 + (D17 (I ,NJ) - P * D18 ( I ,NJ) ) * V(6) 

T(8) = TEN1(I ,NJ) * V(5 ) - D16(I,NJ) * V(6) - V(10) 
T(9 ) = -TEN1(I ,NJ) * V(4) - D15(I,NJ) * V(6) + V(ll) 
T(10)= -PP * D12(I,NJ) * V(3) + PP * D14(I ,NJ) * V(6) 
T (11)= -P * D12(I,NJ) * V(2) - P * D13 ( I ,NJ) * V(6) 
T(12)= -PP * D12CI.NJ) * V( 1) 

GO TO 15 

20 TT1 ( I , L+l , 1 , J) = V(2) 

TT1 ( I ,L+1 ,2 , J) = V(3) 

TT1(I,L+1,3,J) = V(6) 

25 CONTINUE 

DO 30 IJ = 1,12 
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TF1(I,IJ,J) = V(IJ) 

30 CONTINUE 
35 CONTINUE 
HH = HH2 
DO 65 J = 1, 12 
DO 40 K = 1, 12 
V (K) = 0.0 
40 CONTINUE 

V(J) = 1.0 
M = 0 
NJ = 1 

DO 55 L = 1, 50 

45 K = RUNGE (V , T, NJ, M, HH) 

IF(K .NE. 1) GO TO 50 
T(l) = EA2(NJ)*V(12) 

T(2) = V (4) 

T(3) = V(5) 

T (4) = D29 (NJ) * V (8) 4- D291(NJ) * V(9) 

T(5) = D292(NJ) * V(8) - D29(NJ) * V(9) 

T (6) = D21(NJ) * V ( 7 ) 

T(7) = -P * D23(NJ) * V(2) + PP * D24(NJ) * V(3) 

1 + D25 (NJ) * V (4) - D26(NJ) * V(5) 

1 + (D27CNJ) - P * D28(NJ) ) * V(6) 

T(8) = TEN2(NJ) * V(5) - D26(NJ) * V(6) - V(10) 
T(9) = -TEN2(NJ) * V(4) - D25(NJ) * V(6) + V(ll) 

T ( 10)= -PP * D22 (NJ) * V (3) + PP * D24(NJ) * V(6) 
T(ll)= -P * D22(NJ) * V(2) - P * D23 (NJ) * V(6) 
T(12)= -PP * D22 (NJ) * V(l) 


50 

GO TO 45 
TT2 (L+l , 1 , J) = 

V(2) 


TT2(L+1,2,J) = 

V(3) 


TT2 (L+l , 3 , J) = 

V(6) 

55 

CONTINUE 


60 

DO 60 IJ = 1,12 
TF2(IJ, J) = V(IJ) 
CONTINUE 


65 

CONTINUE 



DO 75 I = 1, NPATH 
DO 70 J = 1,3 
DO 70 K = 1,12 
TT1(I,1,J,K) = 0.0 
70 CONTINUE 

TT1(I, 1,1,2) = 1.0 
TT1(I, 1,2,3) = 1.0 
TT 1(1, 1,3, 6) = 1.0 
75 CONTINUE 

DO 80 J = 1,3 
DO 80 K =1,12 
TT2(1, J,K) = 0.0 
80 CONTINUE 

TT2 (1,1,2) = 1.0 
TT2 (1,2,3) = 1.0 
TT2 (1,3,6)= 1.0 
RETURN 
END 
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APPENDIX B 



APPENDIX B 
USER’S INSTRUCTIONS 


This program computes the natural frequencies and the associated mode 
shapes of a nonuniform pretwisted multiple load path rotating rotor blade 
with coupled flapwise bending, chordwise bending and torsional degrees 
of freedom with differential axial displacements in the load paths. 

I DATA CARD 

* NPATH, ISTAGE, I PLOT, NS1, NS2 

* FORMAT(5 I 5) 

NPATH: Number of load paths, maximum = 3. 

ISTAGE: program performs four functions: 

1 - computes the values for frequency determinants only. 

2 - computes natural frequencies and mode shapes. 

3 - computes natural frequencies only. 

4 - computes mode shapes corresponding to given 

natural frequencies. 
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I PLOT: 0 - no mode shape plots. 

1 - plots mode shapes. 

It is preferable to provide data such that the variation of the properties of 
the blade are linear or uniform between any two stations. For a uniform 
blade, it is enough to provide the data at the axis of rotation and the clevis 
for the load paths and clevis and the tip for the blade. 

Note: Load paths - First station should correspond to the axis of rotation 
and the last station should correspond to the clevis. 

Blade - First station should correspond to the clevis and the last 
station should correspond to the tip of the blade. 
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II DATA CARD SET 


* SPAN1, SPAN2, SCH, OMEGA 

* FORMAT(5 E 14.7) 

SPAN1: span of the load paths (inches). 

SPAN2: span of the blade (inches). 

Note: SPAN1 ♦ SPAN2 = Radius of the rotor. 

SCH: semi -chord of the blade (inches). 

OMEGA: rotational speed, (rpm). 

The following data should be provided in the given order. The format is 
(5 E 14.7) unless otherwise specified. 

1. STA1: station locations for load paths, (inches). 

2. STA2: station locations for blade, (inches). 

The following data corresponds to a single load path system: 

3. MASS1 : mass/unit length (Ib-sq.sec/sq.in) . 

4. EIY1: flapwise bending stiffness (Ib-sq.in). 
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5. EIZ1: chordwise bending stiffness (ib-sq.in). 

6. GJ1: torsional stiffness (Ib-sq.in). 

7. El: distance between mass and elastic axis (inches). 

8. BETA1: twist of the load path including the collective (degrees). 

9. KM1S1: mass moment of inertia about the chordwise axis. (Ib-sq.sec). 

10. KM2S1: mass moment of inertia about an axis perpendicular to the 

chordwise axis through the c.g. (ib-sq.sec). 

11. EA1: axial stiffness (lb). 

12. REPEAT THE DATA 3 TO 11 FOR LOAD PATHS 2 AND 3. 

The following data corresponds to the blade: 

13. MASS2: Mass/unit length (Ib-sq. sec/sq. in) . 

14. EIY2: flapwise bending stiffness (Ib-sq.in). 

15. E1Z2: chordwise bending stiffness (Ib-sq.in). 

16. GJ2: torsional stiffness (Ib-sq.in). 

17. E2: distance between mass and elastic axis (inches). 
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18. BETA2: twist of the blade including the collective (degrees). 

19. KM1S2: mass moment of inertia about the chordwise axis. (Ib-sq.sec). 

20. KM2S2: mass moment of inertia about an axis perpendicular to the 

chordwise axis through the c.g. (Ib-sq.sec). 

21. EA2: axial stiffness (lb). 


Ill DATA CARD 


* HI, H, H2 

* FORMAT(3 E 14.7) 

HI: starting frequency (rad/sec). 

H: frequency increment (rad/sec). 

H2: ending frequency (rad/sec). 

CASE 1: ISTAGE = 1 

In this case, only frequency determinants are computed for various 
frequencies. It starts with frequency HI, increments by H until it 
reaches the value H2. 
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/ 

CASE 2: ISTAGE = 2 

In this case, the natural frequencies and mode shapes are computed 
by the frequency scanning technique. Values of frequency 
determinant are scanned starting from the value HI at steps of H 
until the required sign changes are detected or the value H2 is 
reached. So the required number of frequencies or the number of 
frequencies that lie between HI and H2 whichever is less, are 
computed. If two frequencies are closer than increment H, then 
there is a possibility of missing these frequencies. In such cases 
the frequencies are detected by the fact that if any three 
consecutive determinants have the same sign and the absolute value 
of the middle determinant is the smallest of the three then there 
are two frequencies in that range. In this case smaller increments 
are taken to bracket the missing frequencies. However there exists 
a remote possibility of missing these roots also. In such a case a 
closer scanning of the frequency determinants is required. 
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IV DATA CARD 


* FD(I,J) 


* FORMAT(5 E 14.7) 

This card is required if NPATH is greater than 1 and gives the location of 
the j-th load path with respect to the clevis (inches). 


FD(1,J) = h 

y j 

FD(2,J) = h 

z . 

J 

V DATA CARD 

* ITR - number of tension iterations. 


* FORMATCI 5) 


This card is required if NPATH > 1. ITR = 2 will be enough. 
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VI DATA CARD 


* NF 

* FORMAT(l 5) 

This card is required if ISTAGE is not equal to 1. 
NF = number of frequencies, maximum = 10 

VII DATA CARD 

* FREQEN 

* FORMATC5 E 14.7) 

This card is required if ISTAGE = 4. 

FREQEN (J) = natural frequencies. 



VIII DATA CARD 


* BLANK, DOT, STAR 

* FORMAT(3 A 1) 


This card is required if IPLOT = 1. 


BLANK = blank space 
DOT = . 


STAR = * 
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THE NON-ZERO ELEMENTS OF THE STIFFNESS MATRIX ARE: 

ooooooooooooooooooooooooooooooooooooooooooooooooooo 

Kll = -OMEGAS * M * (L + 2.*C9*L2 + C9*C9*L3) 

K22 = 4.*C5*C5*L3 + 12.*C5*C6*L4 + 9 .*C6*C6*L5 

TEM = L - 2.*C5*L3 - 2.*C6*L4 + CS*C5*L5 +2.*C5*C6*L6 +C6*C6*L7 

K22 = C11*K22 - OMEGAS*M*TEM 

K33 = C1*C1*L3 + 2.*C1*C2*L4+ C2*C2*L5 
K33 = K33 * Cll 

K42 = L + C9*L2 + C5*L3 + (2.*C6 + C5*C9)*L4 + 2.*C6*C9*L5 
TEM = C5*L*L + (3.*C6 + 2.*C5*C9)*L3 + 3.*C6*C9*L4 
K42 = (K42 + X J*TEM ) *OMEGAS*M*E*S BETA 

K43 = C1*L3 + (C2 + C1*C9)*L4 + C2*C9*L5 
TEM = C1*L2 + (C2 + C1*C9)*L3 + C2*C9*L4 
K43 = - (K43 + XJ*TEM)*OMEGAS * M * E * CBETA 

K44 = L + 2 . *C9*L2 + C9*C9*L3 

K44 = OMEGAS * (KM2S - KM1S) * CTBETA * K44 

K53 = -C1*L2 + (2*C1*C3 - C2)*L3 

K53 = -(K53 + (C1*C4 + 2.*C2*C3)*L4 + C2*C4*L5)*C11 

K54 = -L2 + (2.*C3 - C9)*L3 + (2.*C3*C9 + C4)*L4 + C4*C9*L5 
TEM = -L + (2.*C3 - C9)*L2 + (2.*C3*C9 + C4)*L3 
TEM = X J* (TEM + C4*C9*L4) 

KS4 = OMEGAS * M * E * CBETA * (TEM + K54) 

K55 = L-4.*C3*L2+2.*(2.*C3*C3-C4)*L3+4.*C3*C4*L4+C4*C4*L5 
K55 = C1I*K55 

K62 = 2.*C5*L2+(3.*C6 + 8.*C5*C7)*L3 + 6.*(2.*C6*C7 + C5*C8)*L4 
K62 = (K62 + 9*C6*C8*L5)*(-C11) 

TEM = L2 + 2.*C7*L3 + (C8 - C5)*L4 - (C6 + 2.*C5*C7)*L5 
- (2.*C6*C7 + C5*C8)*L6 - C6*C8*L7 
K62 ■ K62 - TEM*OMEGAS*M 

K64 = 2.*C7*L3 + 2.*(C7*C9 + C8)*L4 + 2.*C8*C9*L5 

TEM = (L + (C9+4.*C7)*L2 + (4.*C7*C9+3.*C8)*L3 + 3 . *C8*C9*L4)*XJ 

K64 = -OMEGAS*M*E*SBETA*(TEM + K64) 


K66=L+8 .*C7*L2+(16 .*C7*C7+6 .*C8)*L3+24 .*C7*C8*L4+9 .*C8*C8*L5 
TEM = L3+4.*C7*L4+(4.*C7*C7-»-2.*C8)*L5+4.*C7*C8*L6+C8*C8*L7 
K66 = Cll * K66 - OMEGAS*M*TEM 

K71 = - OME GAS*M*C 1 0* ( L2 + C9*L3) 

K77 = -OMEGAS * M * CIO * CIO * L3 

K82 = -C11*(4.*C5*C5*L3 + 12.*C5*C6*L4 + 9.*C6*C6*L5) 

TEM = C5*L3 + C6*L4 - C5*C5*L5 - 2.*C5*C6*L6 - C6*C6*L7 
K82 = K82 - TEM*OMEGAS*M 


C 



K84 =-OMEGAS*M*E*SBETA*(C5*L3 + (2 . *C6+C5*C9 ) *L4 + 2.*C6*C9*L5) 
TEM a XJ*(2.*C5*L2 + (3.*C6+2.*C5*C9)*L3 + 3.*C6*C9*L4) 

K84 = K84 - OMEGAS*M*E*SBETA*TEM 

K86 =2.*C5*L2 + (8.*C5*C7 + 3.*C6)*L3 + 6.*(C5*C8 + 2.*C6*C7)*L4 
+ 9.*C6*C8*L5 

TEH = CS*L4 + (2.*C5*C7 + C6)*L5+(C5*C8 + 2.*C6*C7)*L6 + C6*C8*L7 
K86 = C11*K86 - OMEGAS*M*TEM 

K88 = 4.*C5*C5*L3 + 12.*C5*C6*L4 + 9.*C6*C6*L5 
TEM a C5*C5*L5 +2.*C5*C6*L6 + C6*C6*L7 
K88 a Cll * K88 - OMEGAS * M * TEM 

K93 a -Cll * (C1*C1*L3 + 2.*C1*C2*L4 + C2*C2*L5) 

K94 a C1*L3 + (C1*C9 + C2)*L4 + C2*C9*LS 
TEM = XJ*(C1*L2 + (C1*C9 + C2)*L3 + C2*C9*L4) 

K94 a OMEGAS*M*E*CBETA*(K94 + TEM) 

K95 = -C1*L2 + (2*C1*C3 - C2)*L3 + (C1*C4 + 2.*C2*C3)*L4 
K95 = (K95 + C2*C4*L5)*C11 

K99 = C11*(C1*C1*L3 + 2.*C1*C2*L4 + C2*C2*L5) 

K102 = L2 + C5*L4 + 2.*C6*L5 
TEM = 2.*C5*L3 + 3.*C6*L4 

K102 a OMEGAS*M*E*SBETA*C 10* (K102 + XJ * TEM) 

K103 = -OMEGAS*M*E*CBETA*C10*(C1*L4+C2*L5+XJ*(C1*L3 + C2*L4)) 

K104 = OMEGAS* (KM2S -KM1S ) *C 10*CTBETA* (L2 + C9*L3) 

K105 = 2 . *C3*L4 -L3+C4*L5 + XJ*(2.*C3*L3 - L2 + C4*L4) 

K105 = OMEGAS*M*E*CBETA*C10*K105 

K106 = 2.*(C7*L4 + C8*L5) 

TEM = L2 + 4.*C7*L3 + 3.*C8*L4 

K106 a -OMEGAS*M*E*SBETA*C10*(K106 + XJ*TEM) 

K108 a C5*L4 + 2 .*C6*L5 

TEM a 2 . *C5*L3 + 3.*C6*L4 

K108 = -OMEGAS*M*E*SBETA*C10*(K108 + XJ*TEM) 

K109 = OMEGAS*M*E*CBETA*C 10* (C 1*L4 + C2*L5 + XJ*(C1*L3 + C2*L4)) 
K1010 = OMEGAS * (KM2S - KM1S) * CIO * CIO * CTBETA * L3 


K113 = -C 1 1* (C 1*C3*L3 + (C1*C4 + C2*C3)*L4 + C2*C4*L5) 


K114 = C3*L3 + (C4 + C3*C9)*L4 + C4*C9*L5 
TEM = C3*L2 + (C4 + C3*C9)*L3 + C4*C9*L4 
K114 = 0MEGAS*M*E*CBETA*(K114 + XJ*TEM) 

K11S = Cll*(-C3*L2+(2.*C3*C3-C4)*L3+3 .*C3*C4*L4+C4*C4*L5) 


K119 = C11*(C1*C3*L3 + (C1*C4 + C2*C3)*L4 + C2*C4*L5) 



K1110 = C3*L4 + C4*L5 + XJ*(C3*L3 + C4*L4) 

KlllO = OMEGAS*M*E*CBETA*C 10*K1 1 10 

Kllll = C11*(C3*C3*L3 + 2.*C3*C4*L4 + C4*C4*L5) 

K122 = -C11*(4.*C5*C7*L3 + 6 .*(C6*C7+C5*C8)*L4 + 9.*C6*C8*L5) 
TEM = C7*L3 + C8*L4 - C5*C7*L5 - (C6*C7 + C5*C8)*L6 - C6*C8*L7 
Kill = Kill - OMEGAS*M*TEM 

K124 = C7*L3 + (2.*C8 + C7*C9)*L4 + 2.*C8*C9*L5 

TEM a XJ*(2.*C7*L2 + (3.*C8 + 2.*C7*C9)*L3 + 3.*C8*C9*L4) 

K124 = -0MEGAS*M*E*SBETA*(K124 + TEM) 

K126 =2 .*C7*L2+(8 .*C7*C7+3 .*C8)*L3+18 .*C7*C8*L4+9 .*C8*C8*L5 
TEM = C7*L4 + (2.*C7*C7 + C8)*L5 + 3.*C7*C8*L6 + C8*C8*L7 
K126 = C11*K126 - OMEGAS*M*TEM 

K128 a 4.*C5*C7*L3 + 6.*(C6*C7 + C5*C8)*L4 + 9.*C6*C8*L5 
TEM = C5*C7*L5 + (C6*C7 + C5*C8)*L6 + C6*C8*L7 
K128 a C11*K128 - OMEGAS*M*TEM 

K1210 a C7*L4 + 2.*C8*LS + XJ*(2.*C7*L3 + 3.*C8*L4) 

K1210 a -OMEGAS*M*E*SBETA*C 10*K12 10 

K1212 = 4.*C7*C7*L3 + 12.*C7*C8*L4 + 9.*C8*C8*L5 
TEM = C7*C7*L5 + 2.*C7*C8*L6+ C8*C8*L7 
K1212 a C11*K1212 - OMEGAS*M*TEM 



where 


OMEGAS = tt 2 

; Cl 

= 6/ Z 2 

; C2 

II 

1 

Os 

U> 

; C3 = 2/1 ; 

C4 = 

-3/ £ 2 ; 

C5 = 

= 3/ £ 2 ; 

C6 = 

-2/£ 3 

; C7 = -l/l ; 

C8 = 

1 /a 2 ; 

C9 = 

-l/a ; 

CIO = 

l/£ ; 

Cll = T ; 

L2 = 

l 2 /2 ; 

L3 = 

a 3 /3 ; 

L4 = 

i 4 /4 ; 

L5 = £ 5 /5 ; 

L6 = 

£ 6 /6 ; 

L7 = 

a 1 n 





£ = length of the element ; M = mass per unit length; 
XJ = distance of the element from the axis of rotation; 
E = distance between mass and elastic axes, e; 

SBETA = sing ; CBETA = cos 8 

KM1S = k 2 ; KM2S = k 2 ; CTBETA * cos2g 
m^ 

R 

T = Q 2 mx dx ; = 2,/2 + XJ 

l 

e 


C-4 



